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PREFACE 


This volume of the Annual Status Report describes the results of work performed on the 
development of the boundary element method during the third year of the NASA Hot 
Section Technology program, “3-D Inelastic Analysis Methods for Hot Section Compo- 
nents” (contract NAS3-23697). The goal of the program is to develop computer codes 
which permit more accurate and efficient structural analyses of gas turbine blades, vanes, 
and combustor liners. The program is being conducted under the direction of Dr. C. 
C. Chamis of the NASA-Lewis Research Center. Prime contractor activities at United 
Technologies Corporation are managed by Dr. E. S. Todd. Subcontractor efforts in the 
development of the boundary element method, at the State University of New York at 
Buffalo, are led by Professor P. K. Banerjee. 
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1. INTRODUCTION 


Aircraft powerplant fuel consumption and expenditures for repair /replacement of worn 
or damaged parts make up a significant portion of commercial aviation’s direct operating 
costs. For modern gas turbines, both factors depend heavily on the degree to which 
elevated flowpath temperatures are sustained in the hot section modules of the engine. 
Higher temperatures reduce fuel consumption by raising the basic efficiency of the gas 
generator thermodynamic cycle. At the same time, these elevated temperatures work to 
degrade the durability of structural components (combustor liners, turbine blades and 
vanes, airseals, etc.) that must function adjacent to or within the hot gaspath itself, 
leading in turn to larger maintenance/ material costs. Pursuit of the best compromise 
between performance and durability presents a challenge that will continue to tax the 
ingenuity of advanced gas turbine design analysts for years to come. 

Hot section durability problems appear in a variety of forms, ranging from oxidation/cor- 
rosion, erosion, and distortion (creep deformations to occurrence of fatigue cracking. 
Even modest changes in shape, from erosion or distortion of airfoils for example, can 
lead to measurable performance deterioration that must be accurately predicted dur- 
ing propulsion system design to insure that long-term efficiency guarantees can be met. 
Larger distortions introduce serious problems such as hot spots and profile shifts resulting 
from diversion of cooling air, high vibratory stresses associated with loose turbine blade 
shrouds, difficult disassembly /reassembly of mating parts at overhaul, etc. These prob- 
lems must be considered and efforts made to eliminate their effects during the engine 
design/development process. Initiation and propagation of fatigue cracks represents a 
direct threat to component structural integrity and must be thoroughly understood and 
accurately predicted to insure continued safe and efficient engine operation. 

Accurate prediction of component fatigue lives is strongly dependent on the success with 
which inelastic stress/strain states in the vicinity of holes, fillets, welds, and other dis- 
continuities can be calculated. Stress/strain computations for hot section components 
are made particularly difficult by two factors - the high degree of geometrical irregularity 
which accompanies sophisticated cooling schemes, and complex nonlinear material behav- 
ior associated with high temperature creep/plasticity effects. Since cooling air extraction 
reduces engine cycle efficiency, concerted efforts are made to minimize its use with the 
result that elaborate internal passages and surface ports are employed to selectively bathe 
local regions (airfoil leading edges, louver liner lips, etc.) for which the high temperature 
environment is most severe. These cooling features frequently interrupt load paths and 
introduce complex temperature gradients to the extent that the basic assumptions of one- 
and two-dimensional stress analysis procedures are seriously compromised and the use of 
three-dimensional techniques becomes mandatory. Even in the presence of cooling, com- 
ponent temperature and stress levels remain high relative to the material’s melting point 
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and yield strength values. The combinations of centrifugal, aerodynamic, thermal, and 
other mechanical loadings that typically occur in flight operation then serve to drive the 
underlying material response beyond accepted limits for linear elastic behavior and into 
the regime characterized by inelastic, time-dependent structural deformations. Thus, an 
ability to account for both complexities, three-dimensional and inelastic effects, becomes 
essential to the design of durable hot section components. 

General purpose finite element computer codes containing a variety of three-dimensional 
(brick) elements and inelastic material models have been available for more than a decade. 
Incorporation of such codes into the hot section design process has been severely lim- 
ited by high costs associated with the extensive labor/computer/ time resources required 
to obtain reasonably detailed results. Geometric modeling systems and automated in- 
put/output data processing packages have received first attention from software develop- 
ers in recent years and will soon mature to the point that previous over-riding manpower 
concerns will be alleviated. A prohibitive amount of Central Processing Unit (CPU) time 
is still required for execution of even modest-size three-dimensional inelastic stress anal- 
yses, however, and is chief among the obstacles remaining to be remedied. With today’s 
computers and solution algorithms, models described by a few hundred displacement de- 
grees of freedom commonly consume one to three hours of mainframe CPU time during 
simulation of a single thermomechanical loading cycle. A sequence of many such cycles 
may, of course, be needed to reach the stabilized conditions of interest. Since accurate 
idealizations of components with only a few geometrical discontinuities can easily con- 
tain several thousand degrees of freedom, inelastic analysis of hot section hardware with 
existing codes falls outside the realm of practicality. 

The Inelastic Methods Program addresses the need to develop more efficient and accu- 
rate three-dimensional inelastic structural analysis procedures for gas turbine hot section 
components. A series of new, increasingly rigorous, stand-alone computer codes is be- 
ing created for the comprehensive numerical analysis of combustor liners, turbine blades 
and vanes. Theoretical foundations for the codes feature mechanics of material models, 
special finite element models, and boundary element models. Heavy attention is being 
given to evolution of novel modeling methods that permit non-burdensome yet accurate 
representations of geometrical discontinuities such as cooling holes and coating cracks. 
A selection of constitutive relations has been provided for economical or sophisticated 
description of inelastic material behavior, as desired. Finally, advantages which accrue 
from application of the improved codes to actual components will be demonstrated by 
execution of benchmark analyses for which experimental data exist. 
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2. SUMMARY 


The 3-D Inelastic Analysis Methods program is divided into two 24-month segments: a 
base program, and an option program which was exercised in February 1985. During the 
base program, a series of new computer codes embodying a progression of mathematical 
models (mechanics of materials, special finite element, boundary element) was developed 
for the streamlined analysis of combustor liners, turbine blades and turbine vanes. These 
models address the effects of high temperatures and thermal/ mechanical loadings on the 
local (stress/strain) and global (dynamics, buckling) structural behavior of the three 
selected components. 

The base program (Tasks I and II) dealt with “linear/polynomial” theory in the sense 
that stresses/strains and temperatures in generic modeling regions are linear/polynomial 
functions of the spatial coordinates, and solution increments for load, temperature and/or 
time are extrapolated accordingly from previous information. The computer codes, devel- 
oped during the base program, hereafter referred to as MOMM (Mechanics of Materials 
Model), MHOST (MARC-HOST), and BEST3D (Boundary Element Stress Technology 
Three-Dimensional), are described in the First and Second Annual Status Reports (NASA 
CR-174700 and NASA CR-175060). 

The option program (Tasks IV and V) will extend the models to include higher-order 
representations of deformations and loads in space and time and deal more effectively with 
collections of discontinuities such as cooling holes and coating cracks. Work on Task IV 
has been completed. Continued development of the BEST3D (Boundary Element Stress 
Technology Three Dimensional) code has constituted a very important accomplishment 
of the Task IV effort. The basic BEM (Boundary Element Method) algorithms have 
been extended to allow improved treatment of inhomogeneous elastic material properties, 
thermal loading and plastic deformation. In addition, the time embedded elastodynamic 
algorithm has been expanded to account for nonlinear effects. The results of the boundary 
element method development (Advanced Formulation Development - Task IVC) are given 
in the subsequent sections of this volume. 
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3 . 


BOUNDARY ELEMENT METHOD 


3.1 OVERVIEW 

It is the goal of the advanced formulation development portion of this program to develop 
a computational technique for the solution of linear, nonlinear and transient problems in 
gas turbine engine hot section components. This technique is to be distinct from, and 
complementary to, the Finite Element Method. The existence of such a computational 
method will enhance the ability to calibrate the other codes developed under this contract. 
In addition, it is to be expected that different techniques will prove optimal, in terms 
of efficiency or accuracy, for particular types of component analysis. Since almost all 
general purpose structural analysis computer programs presently available employ the 
displacement finite element method, the new program developed as part of the advanced 
formulation development effort can be expected to extend the ability to perform realistic 
analyses of hot section components. 

During the first year of this program (Task IC), Pratt k Whitney and its subcontractor, 
the State University of New York at Buffalo (SUNY-B), developed a new general purpose 
structural analysis program, BEST3D (Boundary Element Stress Technology), based on 
the use of the Boundary Element Method (BEM). During this work, the boundary ele- 
ment method was implemented for very general three-dimensional geometries, for elastic, 
inelastic and dynamic stress analysis problems. 

In the second year of the program (Task IIC), Pratt & Whitney and SUNY-B continued 
the theoretical and numerical development and the computer implementation of the BEM, 
making very significant advances in a variety of areas, including: 

• Incorporation of substructuring capability in nonlinear boundary element method 
stress analysis. 

• Algorithm development, coding and validation of an embedded time algorithm for 
elastodynamic problems. 

• Initial development of a method for representing the effects of embedded disconti- 
nuities without explicit boundary modeling. 

• Verification of the nonlinear solution capabilities of BEST3D using externally gen- 
erated data. 

The development of the boundary element method during the Base Program is discussed 
in detail in the 1st and 2nd Annual Status Reports (NASA CR-174700 and NASA CR- 
175060). 
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During the first year of the Option Program (Task IVC), the analytical basis of the bound- 
ary element method has been significantly extended and improved capabilities have been 
incorporated in BEST3D. In addition various aspects of the formulation development, 
not yet incorporated in BEST3D, have been tested using a variety of smaller stand-alone 
codes. In particular significant effort has been devoted to the following: 

• Significant improvement of the time embedded elastodynamic analysis, including 
the incorporation of linear time variation of field quantities. 

• The initial development, coding and testing of a boundary element algorithm for 
time-embedded nonlinear dynamic analysis. 

• Significant improvements in the volume integral based plasticity algorithms, in- 
cluding the incorporation of a variable stiffness plasticity algorithm and indirect 
calculation of the most highly singular volume integral terms. 

• Development of a comprehensive formulation for the use of particular integrals to 
represent body force, inhomogeneity and discontinuity effects. 

• Coding and checkout of a (multi-region) BEST3D capability for the determination 
of natural frequencies and mode shapes. 

• Further validation/verification of BEST3D. 

The work carried out in these areas as part of Task IVC is discussed in detail in the 
following sections. 
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3,2 


LITERATURE SURVEY UPDATE 


During the last year there have been two major conferences on boundary element meth- 
ods. One was held in Lake Como, Italy, while the other was a part of the ASME Winter 
Annual meeting held in Miami. It was interesting to note that although nearly 150 papers 
were presented in these two conferences, only a few represented a major advance in the 
state of the art. 

Nishimura described an approximate two-dimensional boundary element formulation for 
the transient coupled problem of flow through a porous medium. Paris and Garrido exam- 
ined the use of nonconforming elements in two-dimensional multiregion contact problems, 
in which they considered opening and closing as well as adhesive and frictional slipping 
and sliding of subregions. Kamiya and Sawaki examined the problem of converting the 
volume integration into an equivalent surface integral using a polynomial expansion of 
the body force field together with the divergence theorem. Li and Brebbia examined the 
use of axisymmetric curvilinear coordinate systems in boundary element analysis in a 
very general manner. All of these papers were presented in the Lake Como conference. 

Ahmad, Banerjee, Wilson, Miller and Snow presented three papers essentially describing 
developments in inelastic and dynamic analysis as incorporated in BEST3D. Rizzo and 
Shippy described developments in periodic dynamic analysis, similar to methods incor- 
porated in BEST3D. Mukherjee and Chandra recently presented further developments 
in large deformation analysis of axisymmetric solids by the boundary element method. 
Brebbia reported a novel boundary element method formulation for transient and periodic 
two-dimensional problems using a real variable method which utilizes the fundamental 
solutions of the static elastic problem. Annigeri, Gerstle and Ingraffea, Polch, et. ai 
examined a variety of problems involving the modeling of crack tip singularities. Kline in 
an excellent paper examined the use of parallel processing and solution of block banded 
system equations in a substructured boundary element analysis. All of these papers were 
presented in the ASME Winter Annual Meeting (Cruse, T. A., A. Pifko and H. Armen, 
Advanced Topics in Boundary Element Analysis , AMD Vol. 72, ASME, N.Y., 1985). 
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3.3 


FORMULATION DEVELOPMENT 


3.3.1 SUMMARY 

Important advances have been made during Task IVC in extending the BEM formulation 
for three-dimensional stress analysis of gas turbine engine structures. The most signif- 
icant advances have been in the areas of dynamic analysis, where an improved linear 
formulation and an initial nonlinear formulation were developed and implemented, in the 
exploitation of particular integrals of the equilibrium equations to represent a wide variety 
of volume based phenomena and in the three-dimensional developmentof a variable stiff- 
ness plasticity algorithm, incorporating the constitutive relationship within the boundary 
element system matrix. The extended or newly developed formulations are discussed in 
the subsections below. The basic BEM formulation is only very briefly reviewed, as full 
details are available in the First and Second Annual Status Reports (NASA CR-174700 
and NASA CR-175060). 

3.3.2 QUASI-STATIC ANALYSIS 

3.3.2. 1 REVIEW By making use of the reciprocal work theorem, the governing 
differential equations for a three-dimensional (homogeneous) structure under combined 
thermal, mechanical and body force loadings can be converted to an integral equation 
written on the surface of the structure. This integral equation is: 

cytii = / (GijU - FijUi) dS + [ ( Gufi + (3WjT) dV (3.1) 

J s Jv 

where T = temperature, Wj = TikjSik, (3 = coefficient of thermal expansion, T^j = the 
stress, cr,fc, due to a point force system, ej, and G^, T,*, and are defined in the First 
Annual Status Report (NASA CR-174700). The equation 

&ij = Jj^Dijktk — SijkUk) dS + J^(T x j k t k + MijT)dV (3.2) 

allows calculation of stresses at any interior point where they are required. A similar 
equation for interior displacements can be obtained by setting = 8 XJ in (3.1). 

In a purely elastic problem BEM stress analysis can be carried out entirely on the bound- 
ary of the structure. Once a physically reasonable set of boundary conditions has been 
prescribed, (3.1) can, in principle, be solved for all of the remaining boundary displace- 
ments and tractions. 

It is generally impossible to solve (3.1) exactly for real structures and loading conditions. 
Suitable approximations of the boundary geometry, displacements and tractions must be 
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used in order to reduce (3.1) to a system of algebraic equations. The present version 
of BEST3D models boundary geometry and boundary values of field quantities using 
linear and/or quadratic isoparametric shape functions. The surface integrals in (3.1) 
are then evaluated numerically using product Gaussian quadrature rules. The numerical 
implementation of the BEM is discussed in detail in textbooks by Banerjee and Butterfield 
and by Brebbia, as well as in the First and Second Annual Status Reports (NASA CR- 
174700 and CR-175060). 

In the case of inelastic analysis, the volume integrals in (3.1) cannot be calculated a priori , 
since they require knowledge of inelastic strain, which is itself a part of the solution. In 
this case equations (3.1), (3.2) and the inelastic material model can be regarded as a 
coupled system of nonlinear equations. In the numerical implementation of the BEM 
(3.2) is used to calculate the stresses at interior points, and the nonlinear material model 
is then used to evaluate inelastic strain. Since the volume integrals of inelastic strain 
vanish except in regions of nonlinear material response, approximations of geometry and 
field quantities are required only where nonlinearity is expected. In BEST3D volume 
models utilize quadratic isoparametric shape functions to model strain variation within 
volume cells. Three constitutive models are presently included in BEST3D: 

1. Von Mises model with isotropic variable hardening 

2. Two-surface model for cyclic plasticity 

3. Combined plasticity and creep model of Walker. 

These constitutive models are discussed in detail in the First and Second Annual Status 
Reports (NASA CR-174700 and NASA CR-175060). 

The remainder of Section 3.3.2 discusses significant new developments carried out during 
Task IVC. 


3. 3. 2. 2 VARIABLE STIFFNESS PLASTICITY A new approach is outlined for 
BEM formulations for elasto-plasticity, which exploits certain features of the constitutive 
relationships involved. The unknown nonlinear terms in the interior are now defined 
as scalar variables. A new direct numerical solution scheme comparable to the variable 
stiffness method used in finite element analysis has been developed and applied to a 
number of standard plasticity problems. 

For a standard elasto-plastic flow problem the evolution of plastic flow is governed by: 


F(eij,h) = 



dF 


'da. 


‘j 


(3.3) 

(3.4) 
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These equations together with the consistency relations (i.e., the stress point must remain 
on a newly developing yield surface characterized by a change in the hardening parameter 
h) leads to an expression for the unknown plastic flow factor A as: 


A = L' <7, 




(3.5) 


where 


, _ 1 OF 
ij ~ H dcrij 


*--<#■ 


d_F dh dF 
dh d *Ln d<r ™ 


It should be noted that L depends upon the current state variable, not on the incre- 
mental quantities. 


However, the relationship given by (3.5) does not exist for ideal plasticity, as H vanishes 
for zero hardening. This can be avoided by reformulating the above expression in terms 
of strain increments: 

A = LUa (3.6) 

where 


V.. 




H' 


l_ dF 
H' d °ij 

dF n dF df 

da kl Uklmn damn [ de p kl + dh de P kl > da kl 


where Diju is the elastic constitutive tensor. It is evident that H l does not vanish for 
zero hardening (ideal plasticity). 

The basic boundary element formulation for an inelastic body undergoing infinitesimal 
strain is given by: 


ciMO = JjGifatM *) - £;(*. 0*(*)J ^ <7° dn (3.7) 

The stress rates at an interior point £ are obtained from equation (3.7) via the strain- 
displacement relations and the constitutive relationships (d tJ = Z)* ) as 

■>«(«) = ( 3 - 8 ) 

+ /„ **(.) da + r, k ^(0 


where the kernel functions have been defined in the Second Annual Status Report (NASA 
CR-175060). 
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In equation (3.8) the volume integral must be evaluated in the sense of (fi — D ) with the 
limit taken as D — ► 0, where D is a spherical exclusion of small arbitrary radians with £ 
as its center. The term J " is the jump term derived from the analytical treatment of the 
integral over D. It is of considerable interest to note that the value of J * is independent of 
the size of the exclusion D, provided the initial stress distribution is locally homogeneous, 
i.e. uniform over its volume. 

The evaluation of strains and stresses at boundary points can be accomplished by consid- 
ering the equilibrium of the boundary segment and utilizing constitutive and kinematic 
equations. The stresses and global derivatives of the displacements which lead to strains 
at a point £ can be obtained from the following set of coupled equations: 

+ + <»(£))) = -*°(0 

= U(t) (3-9) 

d&duiH 1 _ diiitt) 

dm dik diji 

where rji are a set of local axes at the field point £. 

All the above nonlinear formulations include initial stresses in the governing equations 
which are not known a priori and, therefore, are solved by using iterative procedures. 
A non-iterative direct solution procedure is made feasible in this work by reducing the 
number of unknowns in the governing equations by utilizing certain features of the incre- 
mental theory of plasticity expressed by equations (3.3) to (3.6). The initial stresses <r°. 
appearing in equations (3.7) to (3.9) can be expressed in the context of an elastoplastic 
deformation as: 

= KijX (3.10) 

where K u = D’ ju ^. 

Substituting (3.5) and (3.10) in equations (3.7) and (3.8) we can obtain: 

dMO = f r [Gii(*)U(x) ~ Fij(x,()ui(x)]dT (3.11) 

and 

ka = {)*.(*)] <ir 

+£'»(«) jf B ’ P . ,(*,{) Ki p {x) \(x) da (3.12) 

+f'„ •>:„* *»«) m 

Equations (3.11) and (3.12) can be solved simultaneously to evaluate the unknown values 
of displacements, traction rates and the scalar variable A. 
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The equations for the boundary nodes (3.9) are similarly transformed to express them in 
terms of the scalar variable A using equations (3.5 ) and (3.10). 


3. 3. 2. 3 VOLUME INTEGRATION ALGORITHM Assuming that the bound- 

ary and a representative area of the volume (which is likely to yield) is divided into 
N isoparametric boundary elements and M isoparametric volume cells respectively, we 
can express the coordinates as well as the displacements, tractions and the scalar flow 
parameter A in terms of their nodal values as: 


Xi = N a {vi)x a i 
m = N a (v)u* 

U = N a ( V )i~ 

over surface patches, and 

x, = M 

<x° = M (3 (7j)a J3 . 

Xi = M 0 ( v )K t] >? 


over volume cells, where a and / 3 represent, respectively, the nodal points of boundary 
elements and volume cells. 


The boundary integral equations for displacement (3.11) can then be represented as: 


C t ,Ui(x P Q ) 


!){/ G ii [ X n (T,),x P 0 }N a (T 1 )t™dY n 

n = 1 JVrx 

- I f'aW'i'n),* o] n ‘ x ( t i ) C dTn ) 

J r n 



Bi„[x m ( V ),x P 0 ] M 0 (t]) Ki p \ 0n dQ, 


(3.13) 


The evaluation of the surface integrals involved has been extensively documented in 
the published literature, notably by Lachat and Watson, who first developed the use of 
isoparametric surface elements for elastic analysis. A more efficient technique developed 
in this work has been outlined in the previous annual report (NASA CR-175060). 

The volume integral is transformed and the quadrature rule is applied as: 



(3.14) 


= JJJ BiAzil)’ £l M 0(v) J (v) di) I dr )2 drj 3 
= £ 't'tw a w b w c B ljk [x,T f abc ,t] M 0 ( V abc ) J( V abc ) 

a= 1 b— 1 c=l 

where the volume element and Jacobian are defined by: 

dClm = J(vi . m , Vi ) dvi dr )2 d-q 3 
( 3x * \ 

= i,j = 1,3 

The non-singular integration is then carried out in a straightforward manner using ordi- 
nary Gaussian integration formulae, with appropriate error control via cell subdivision. 

The singular integration technique is selected by carefully studying the behavior of the 
kernel functions involved. For the volume integral in (3.13 ) which is singular of the order 
1 /r 2 , the adopted procedure of integration is as follows: 

1. The curvilinear cell is transformed to a unit cell by using the isoparametric trans- 
formation. 

2. Depending upon the location of the singular point, the cell is subdivided into three 
sub-cells, when the singularity is at a corner node, or four sub-cells, when the 
singularity is on the cell side. 

3. Each sub-cell is then mapped into a unit cube. 

Since the original integrand is singular of the order l/r“, it becomes bounded in the 
transformed domain and, therefore, it can be evaluated accurately by selecting a suitable 
order of integration rule. 

The singular integral involving the function Bi p jk is strongly singular (of order 1/r 3 ) and 
the integrand is made bounded by excluding a neighborhood of the source point, in the 
form of a sphere, and mapping the remaining curvilinear cell to a unit cell for integration 
using the quadrature formula. This cell is further sub-divided in a manner which reflects 
the order of the integrand involved. 

A much more satisfactory procedure can be used to evaluate the singular contributions in 
the interior stress calculations, if an entire subregion is filled with cells. If the subregion 
is unrestrained, the surface tractions will be zero and the displacements can be easily 
calculated for the case of uniform initial stress distribution. The singular terms can be 
then easily calculated by summing up the contributions such that the stresses are zero 
at every nodal point. 


12 


In the present version of BEST3D both above mentioned procedures have been imple- 
mented. 


3. 3. 2. 4 ANISOTROPIC ELASTICITY Anisotropic materials are frequently 
called for in the design of the hot section of a gas turbine engine. The material anisotropy 
typically arises either from a controlled casting operation (as with directionally solidified 
or single crystal airfoils) or from processing operations applied to an originally isotropic 
material. The anisotropic material classes presently of major interest in the analysis of 
gas turbine engine structures are transversely isotropic (directionally solidified castings) 
and cubic (single crystal castings). In addition, local nonlinear deformation may induce 
anisotropy in an otherwise isotropic material. 

As pointed out previously (Ref. 2) the application of the boundary element method to 
the elastic analysis of anisotropic materials requires, as a minimum: 

• The point load solution for the (generally oriented) anisotropic material must be 
calculable. 

• The surface stress calculation must be appropriately modified. 

In addition, application to meaningful problems requires the development of a particular 
integral for centrifugal loads. Further, if nonlinear analysis of such materials is to be 
considered, it must also be possible to calculate the higher order kernels required for the 
evaluation of interior stress and for the volume integration of initial stress. 

The BEST3D surface stress routines were originally written to handle general anisotropic 
materials, and the required particular integral for centrifugal loading was developed dur- 
ing Task IIC. The major challenge in the efficient application of boundary element analysis 
to anisotropic materials is thus the development of suitable techniques for the calculation 
of the required kernel function values. 


3. 3. 2.4. 1 Analytical Kernel Representation The point load function (displace- 
ment kernel) for a general anisotropic material can be represented as a line integral, 

This equation defines the displacement in the i-direction at a point x due to the appli- 
cation of a point load in the ^-direction at the point y. The line integral is taken on 
the circle on the unit sphere which is perpendicular to the vector x—y. The integrand, 
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defined by, 


K,M) = extern (3.16) 

is the inverse of a matrix, whose entries are quadratic forms with coefficients which are 
the elastic constants of the (Cijkm) of the material. 

The ideal solution to the calculation of the point load function, and the required deriva- 
tives, is the reduction of the line integral to a closed form, elementary expression. In the 
case of an isotropic material, \Kij\ is constant over the unit sphere. As a result the line 
integral simplifies considerably, leading finally to the familiar Kelvin point load solution. 

In the case of anisotropy, the evaluation of (3.15) requires further consideration. Using 
Fourier transform techniques (Ref. 3), it is possible to develop an alternative represen- 
tation of the point load solution as an integral in the complex plane, 

Uij{* -y)= l ~ dz = Uij(r) (3.17) 

where 2 is derived from £ by a change of variables to spherical polar coordinates. The 
denominator T{z) is a 6th order polynomial equation in 2. Complex variable theory allows 
the representation of the point load solution in terms of the residues of the integrand at 
the zeroes of T(z) (the poles of the integrand). 

In order for this representation to be useful in boundary element analysis it is necessary to 
solve the equation T{z) = 0 explicitly, in terms of the elastic constants of the material, for 
an arbitrary orientation of the vector x — y. This can be done, trivially, for an isotropic 
material. For the important case of a directionally solidified material the problem is 
much more complex, but still solvable. In that case T[z) factors into the product of 
three quadratic expressions. The roots of each factor can therefore be calculated in terms 
of radicals, allowing the use of the residue theorem to produce a closed form solution in 
terms of elementary functions. A particularly useful form of this solution was developed 
in Ref. 4, and was, after considerable simplification, incorporated in BEST3D during 
Task IIC. The same authors have extended this work to produce a closed form, but very 
complicated, point load solution for two-phase transversely isotropic materials (Ref. 5). 

3. 3. 2.4. 2 Single Crystal Materials The other material class of primary interest 
in hot section components consists of cubic crystals (single crystal materials). These 
materials are defined by three material constants (Cn, C 12 , C44). Despite this seemingly 
simple structure , T{z) can be solved using radicals only for certain very special cases: 

• Cu — C12 — 2C44 = 0 (isotropic) 

• C12 = -C44 
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• x — y is aligned with a material axis. 

The last of these special cases is not useful, because the evaluation of the point load solu- 
tion is required for general orientations of x — y. Single crystal turbine engine materials 
satisfy neither of the first two conditions. 

In order to apply boundary element analysis to single crystal materials it is, therefore, nec- 
essary to develop techniques for the efficient and accurate numerical approximation of the 
point load solution and its derivatives. There are three major categories of approximate 
representations which can be considered: series representations, tabulation/interpolation 
of kernel function values and matrix techniques. These techniques will be reviewed in 
the context of their suitability for boundary element analysis. 

3. 3. 2. 4. 3 Series Representation In a series representation each term of the point 
load solution is represented as an infinite series of relatively simple functions, frequently 
the powers of an expansion parameter. A frequently used parameter for single crystal 
materials is d = C xx — C X 2 — 2 C 44 , the ’deviation’ from isotropy. The coefficients in the 
expansion are chosen to give a good fit to the actual point load solution using only a few 
terms of the series representation. Such techniques are discussed in detail in Reference 3. 

Three major disadvantages make these techniques unsuitable for use in boundary element 
analysis. First, accuracy is difficult to evaluate or control. In general the success of the 
approximation depends on how near the single crystal material is to an isotropic material. 
As Figure 3.1 indicates, aerospace single crystal alloys are relatively far removed from 
isotropic behavior. Second, the algebraic process required to evaluate the coefficients 
in the approximation is very complex. Finally, the derivatives of such representations 
typically lose further accuracy with every differentiation. Depending on the application, 
boundary element analysis requires from one to three derivatives of the basic point load 
solution. It would not be possible to maintain required accuracy using derivatives of a 
series representation. 


3. 3. 2 . 4. 4 Interpolation in Tabulated Data A second approach to the approxi- 
mation of the single crystal point load function is the use of tabulated data together 
with suitable interpolation routines. This approach is described in some detail in Ref. 6. 
Briefly, it is possible to numerically evaluate the line integral in (3.15) to any required 
degree of accuracy. Further, the line integral depends only on the orientation of x — y. 
It is thus possible to tabulate the values of the line integral on the surface of the unit 
sphere, and to interpolate in the table to obtain the orientation dependence of interme- 
diate points. Further, by differentiation of (3.15) and repeated use of the chain rule, it 
is possible to reduce any required derivative of the point load function to related, non- 
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Figure 3.1: Relationship of Single Crystal Material to Isotropy 


singular orientation dependent line integrals. These integrals can also be tabulated, with 
the same degree of accuracy as the point load function. 

The method outlined has been implemented in a three-dimensional boundary element 
code, with generally excellent results. Analysis is significantly more costly than for an 
isotropic material because of the large amount of computation required for the inter- 
polations. The penalty in computing time is a factor of approximately 1.9 compared 
to approximately 1.3 for the closed form directionally solidified kernel implemented in 
BEST3D. 

In addition to the question of computing time the tabulation/interpolation technique has 
major disadvantages for use in a general purpose code, such as BEST3D. First, a great 
deal of data is involved. Using a 17x33 grid to cover the unit sphere, the point load table 
requires almost 3400 (double precision) words of storage. The first derivatives (required 
even for the boundary solution) require an additional 10,000 words. Second derivatives 
(required for interior stress calculation) require 30,300 words. Third derivatives (required 
for volume integral based plasticity) would require 90,000 more words. Second, the tables 
(produced and stored by a separate program) are based on a specific set of material con- 
stants. If the material properties are changed, then the entire table must be regenerated. 
Finally the questio of temperature dependent material properties is extremely difficult. 
Data storage requirements are, of course, multiplied by the number of temperature values 
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in the data base. More importantly, there is no reliable method for interpolating in the 
data base to obtain kernel function values at intermediate temperatures. 

3. 3.2,4. 5 Matrix Formulation An alternative formulation can be given in terms 
of matrix algebra. This alternative formulation was initially developed in the field of 
solid state physics, in a context much more general than the simple construction of static 
point load solutions. The general technique can be applied to the complete equilibrium 
equation, including the inertia term. Further, more general singular distributions (for 
example, line loads and dislocations) can be considered, as can the presence of initial 
plastic strain. The possible application of the technique for the calculation of the single 
crystal point load solution is described below. The notation employed is generally that 
of K. Malen (Ref. 7). 

In the static case, the equilibrium equation takes the form, 

dijj ~ ~ft (3.1 8) 


with the stress-strain relation, 
and 


&ij = Cijkifiki 


U k,i — fiki 


Taking Fourier transforms and substituting for /3^ in terms of £/*, the following equations 
are obtained: 


%}cj 


= -pwUi-fi 

= ikiCijkiUk 


(3-19-) 


To simplify the calculation of the point load solution a new coordinate system is intro- 
duced. The vector r is a unit vector in the direction from the load point to the field 
point. The unit vectors m and n are then chosen to complete a right handed orthogonal 
coordinate system. The vector from the load point to an arbitrary point in space can 
then be represented as: 

k = k(ccm + fin + 7 r) 


with, 


k = |fc| 


Adopting the convention that 
defining 


(**)ifc — kjCjikih 


p = f3/a 
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q = 7 /a 

and taking E to be the 6 x 6 identity matrix, the equilibrium equations can be written in 
matrix form as 

(N -pE + qAX = k~ 2 a~ 2 Q (3.20) 

where the vectors ( and Q are defined as 

Q = Ui Qi = 0 *=1,3 

C* = ri j&tj Qi — ifi z — 4,6 


and _ 

N _ — (nn)" nm —i(nn) 

+ i(mn)(nn)~ l (nm) — (mn)(nn ) _1 

A _ [ -(nn) _ 1 (nr) _ 0 

— t[(mr) + (rm) + p(rn) + g(rr)] + z(mra)(nn) (nr) 0 


In the case of a point load, fi = FiS(r), so that, 

h = * 

and the vector Q is independent of k, simplifying to 

Qi = 0 i = l,3 
z = 4 , 6 


leading to the equation 


aP,q) = (N- P E + qA)- l Q 


(3.21) 


Fourier inversion leads to an expression for (, in physical space in terms of the eigenvalues 
of (N-pE). 

f/ r \ - Jl-Y' i \ Adj(N — pE)Q ] ('3 22') 

The first three components of this vector are the Cartesian components of displacement 
at r due to the application of the point load at the origin. The plus sign in this equation is 
used when G(p v ) > 0 and the minus sign when $s(p v ) < 0 . While there is a considerable 
amount of calculation required for the evaluation of this expression, all of it involves 
straightforward matrix operations. Even the extraction of the eigenvalues of the 6 x 6 
matrix ( N — pE ) can be performed using standard software packages. Further, in the case 
of the single crystal material it will be possible to considerably simplify the calculations 
involved in this evaluation. Evaluation of higher order kernels will be carried out by 
differentiation of the matrix form of the displacement kernel. Although additional terms 
are generated , the basic eigenvalue problem remains the same for the higher order kernels. 

The utilization of the matrix formulation for evaluation of anisotropic point load solutions 
has several major advantages: 
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1. The need for the storage and manipulation of large, pre-calculated kernel function 
tables is eliminated. 

2. Kernel function evaluations can always be done using correct material properties, 
eliminating the need to define a temperature interpolation at the kernel function 
level. 

3. Accuracy will be easier to control, since no numerical integration is required in the 
kernel function evaluation. The only potentially sensitive operation is the solution 
of the 6x6 eigenvalue problem, but the currently available algorithms allow this to 
be done with very high accuracy, without user or programmer intervention. 

It is anticipated that the calculations involved can be streamlined sufficiently to allow 
the exact calculation of kernel function values at all required points. In the event that 
this cannot be done, it is planned to calculate the exact values, for a given source point 
location, at a suitable set of points (normally the geometric nodes) on a single surface 
element. Interpolation using the normal isoparametric shape functions will then be used 
to provide a good approximation to the kernel function values at individual integration 
points. 


3. 3. 2. 5 APPLICATION OF PARTICULAR INTEGRALS The early devel- 
opment of the boundary element method, especially in three dimensions, concentrated 
on the static elastic analysis of homogeneous materials. Recent developments, including 
work done in the present program, have demonstrated that the method is applicable to 
both anisotropic elastic analysis and inelastic analysis. In particular the present version 
of BEST3D includes provision for the analysis of directionally solidified materials and 
for inelastic analysis. It is planned to incorporate the kernel functions for single crystal 
materials during the final year of the program, as outlined in 3. 3. 2. 4. 

The major remaining technical obstacle to the incorporation of the three dimensional 
boundary element method in the aerospace design process is the difficulty encountered in 
efficiently handling continuously inhomogeneous (generally temperature dependent) elas- 
tic material properties, and, to a lesser extent, thermal stress. Multi-material problems 
can be effectively handled by substructuring, or (potentially) the use of special kernel 
functions. Situations with continuously variable elastic properties and/or time depen- 
dent temperature distributions have required the use of a subsidiary volume discretiza- 
tion and volume integration. These procedures are discussed in the first two Annual 
Reports (NASA CR-174700 and NASA CR-175060). A two-dimensional application to 
elastic inhomogeneity is presented in Ref. 8. 

While this technique is capable, in principle, of providing the solution to very general 
three-dimensional problems it is not suitable for routine application in an aerospace design 
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environment, for two major reasons: 


1. The creation of the volume discretization, while less demanding than in finite ele- 
ment analysis, still imposes a significant additional burden on the analyst. 

2. The extensive volume integration involved is extremely expensive. Further, despite 
significant advances made during Task IVC (3. 3. 2. 3) it remains very difficult to 
guarantee sufficient accuracy in the volume integrals, particularly when thermal 
loading is important. 

Plasticity analysis is also affected by the costs associated with volume integration and 
volume discretization, although the cost of finite element plasticity analysis makes the 
expense associated with volume discretization in the boundary element method appear 
less significant. In addition, in many plasticity analyses the region of plastic behavior is 
quite constrained, helping to control both the size of the volume mesh and the associated 
costs. 

The method developed during Task IVC is based on the well known fact that the solution 
to an inhomogeneous system of differential equations can be represented as the sum of a 
complementary solution satisfying the homogeneous equation and a particular solution 
satisfying the inhomogeneous equation. In the case of the stress analysis of a solid the 
equilibrium equations can be written as: 

+ (A + - Af“ 4i . + 2/zcV. (3.23) 

or in operator form as: 

= ' 4 , + 2 * 4 ,- ( 3 ' 24 ) 

where A and /x are reference moduli (constant throughout a substructure) and e 0 .. is 
the initial strain, generally nonzero throughout the part, due to thermal, inelastic and 
inhomogeneity effects. The initial strain is written as: 

«£(*) = £*«(*)*«(*) - + &H&T(x)a(x) (3.25) 

where if, H° and #(= H — H°) are the reference, local and deviation compliance tensors, 
respectively. If a particular solution to (3.24) can be found, then the physical solution 
can be represented as: 


u 


physical complementary 


— U 


+ vP' 


>ar ticular 


and a boundary integral equation can be derived for the complementary solution. Ex- 
pressing the displacements and tractions of the complementary solution in terms of the 
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physical solution and the (known) particular solution then allows the development of a 
boundary equation for the physical solution: 

CijUj - JjGijti - FijUi] dS = CijV?. - J [Gijf. - Fijuf ] dS (3.26) 

This is then reduced to algebraic form in the usual manner. A similar technique is already 
employed in BEST3D for the solution of problems with centrifugal loads. 

Application of this method to the solution of problems of inhomogeneity and thermal 
stress will depend on the development of suitable particular solutions to be used for the 
representation of the initial strain fields involved. 

It can be shown (Ref. 9) that (assuming appropriate continuity and differentiability 
conditions), if a potential function g is chosen so that: 

0 — 

^ — Qijykkmm 

then a Galerkin vector F, defined by: 

(1 - u)F t = \g k k,, + 2 figijj 

will lead to a displacement field which is a particular solution of (3.24). Any solution 
of the inhomogeneous biharmonic equation thus becomes a possible starting point for a 
particular solution of (3.24). There are two general classes of suitable functions: 

1. Domain functions, chosen to allow representation of a suitable variation of initial 
strain throughout an entire substructure or part. 

2. Point based functions, which take on their largest value at a particular boundary 
or interior point and decay as a function of the distance from that point. 

Domain functions are particularly suitable for the representation of thermal strain, whose 
variation is usually known a priori, while point based functions appear more suitable for 
the representation of initial strain due to inhomogeneity. If necessary, both types of 
representations can be used in a single analysis. 

Both forms of initial strain approximation discussed above are continuous and nonzero 
everywhere in three-dimensional space. An additional class of representations exists, in 
which the nonzero initial strain distribution is restricted to a small volume of space. This 
class is particularly useful in the representation of holes, without explicit modelling. It 
will be discussed further in 3. 3. 2. 6. 

The remainder of this section outlines one possible algorithm for the treatment of the 
combined thermal and mechanical loading of a structure with temperature dependent 
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elastic properties. In this case the initial strain is: 

e“. = H ljk ,<7ki + SijaAT (3.27) 


Suitable functions are chosen for the approximation of the thermal strain (typically do- 
main functions M ni with weights V*")* and the inhomogeneity strain (typically point 
based functions M mi with weights ). The approximation of initial strain is written as: 

«?,• = £ M n i>” + £ (3.28) 

n— 1 m— I 

It should be noted that the sums in the above equations may be either over sets of points 
(surface and interior) or over a set of terms not specifically associated with geometric 
locations in the structure. 


The initial strain approximations lead, through appropriate differentiation and use of 
the stress-strain relationship, to particular solutions for the thermal and inhomogeneity 
initial strains. It is important to note that the particular solutions are linear in ip n - and 
0™, respectively . 


If the physical solution is represented as: 


! 

u = u - h t/ 


then the usual methods allow the derivation of the algebraic form of the boundary integral 
equation for the complementary solution: 


Ax c = By 


where z is the vector of unknowns and y_ the vector of known displacements and tractions. 
This can be rewritten in terms of the physical unknowns as: 


Ax = By + AQ<f)_ 4- AQiy 


(3.29) 


A similar development allows representation of the stress, at interior or boundary points, 
as: 

o_ — Cx_ + Dy + C R(j> + C 

The vector of initial strains due to inhomogeneity can then be written as: 

e° = Ha = H(Cx + Dy c + C R± + CR±) 


and since, 


$ — Pt 
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the equation 


(3.30) 


± = PH(Cx + Dy c + CRl + CRf) 

is obtained. 

These equations can be written in partitioned matrix form as: 

’ A AQ 
PHC PHCR 


X 


rx 



r 2 


In this system the matrices A and C are the usual static boundary element matrices. All 
of the other matrices involved are related either to the stress-strain relationship or to the 
initial strain approximations, and require no integration, only function evaluations, for 
their construction. 

The linear equation system can either be solved directly or the weights, <£, can be elimi- 
nated as a first step, leading to variable stiffness formulation. 

An entirely similar formulation can be developed for nonlinear analysis, including thermal 
effects. Initial versions of such a nonlinear algorithm have been developed and are being 
evaluated and refined in BEST3D. 

In summary, the use of particular integrals now offers the realistic possibility of analyzing 
structures with thermal loading, inhomogeneity and nonlinear response. Input definition 
is significantly simplified, since no interior connectivity is required. Since volume integra- 
tion is no longer required, the highly singular 4th rank kernels are completely avoided, 
as is the (highly singular) volume integration of the third rank kernels. The only direct 
attention required to the interior of the structure is the possible placement of points for 
use in the approximation of initial strain. 


3. 3. 2. 6 REPRESENTATION OF HOLES/CAVITIES Turbine engine hot sec- 

tion components frequently contain arrays of holes, usually to allow the passage of cooling 
air. On occasion arrays of holes are also used to reduce component weight. In addition, 
cracks may be found in such components, frequently originating from cooling holes or 
other structural details. 

Clearly, one method of solution for such problems is the modelling of the hole or crack 
surface(s) using boundary elements. This method has been shown to accurately define 
stress concentration and fracture mechanics parameters in complex three-dimensional 
structures. Good agreement has been demonstrated between nonlinear experimental and 
boundary element results in notches (Ref. 2). This direct modelling method has the 
advantage of providing very detailed and accurate results near the hole or crack. If only 
a single such discontinuity is to be studied, the additional analysis cost is usually not 
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excessive. If, however, an array of discontinuities is present in a part, direct modelling 
will no longer be feasible. 

Three alternative methods are under active consideration for incorporation in BEST3D 
for the representation of holes and/or cracks without detailed boundary element mod- 
elling, a hybrid surface integral/boundary element method, the development of approxi- 
mate kernel functions and the use of exact or approximate solutions for inclusions. 


3. 3. 2. 6.1 Surface Integral/Boundary Element Hybrid The surface integral me- 
thod (Ref. 10 and 11) allows the modelling of cracks in an unbounded solid using a 
distribution of point force dipoles on the crack surface. Dipole strengths are assigned at 
suitably chosen points on the crack surface (not required to be planar). The resulting 
integral equations are numerically solved at a set of collocation points on the crack surface. 
While points must be defined on the crack surface, no connectivity definition is required. 
The upper and lower crack surfaces are not modelled separately. The surface integral 
method has been successfully coupled to the finite element method for the solution of 
problems involving cracks in finite bodies (Ref. 12). The coupling involves a finite 
element solution for the uncracked body and a surface integral solution for the actual crack 
geometry in an unbounded solid. The two solutions are coupled by enforcing displacement 
compatibility and equilibrium on the external boundary and the crack surface. 

A similar coupling between the surface integral and boundary element methods is clearly 
feasible. In this case the coupling will be made by calculating crack surface tractions 
using the boundary element interior stress identity. 


3. 3. 2. 6. 2 Approximate Kernel Functions A second approach involves the con- 

struction of approximate kernel functions which include the presence of one (or more) 
holes or cracks. These approximate kernel functions can be generated by the following 
process: 

1. A standard boundary element model is built for the desired cavity shape. The 
system matrix corresponding to the exterior problem is built and decomposed. 

2. The tractions on the cavity surface due to an arbitrarily located point force are 
calculated. The negative of these tractions is applied to the cavity surface, and 
the decomposed system matrix is used to solve for the displacements on the cavity 
surface. 

3. The boundary element interior displacement and stress identities are used to cal- 
culate the point load solution at an arbitrary point in the exterior of the cavity. 
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The considerable geometric symmetry which exists for many common hole shapes (cylin- 
der, sphere, ellipsoid) can greatly simplify the calculations. The greatest potential dis- 
advantages of this approach are its numerical complexity and the difficulty of controlling 
accuracy when the source point is near the hole. The method would become quite cum- 
bersome if multiple holes were required in a single subregion. This method does, however, 
share with the surface integral method the ability to account directly for the variation of 
stress in the uncracked part within the region occupied by a hole. 

3. 3. 2. 6. 3 Inclusion Solutions A final method, particularly suited to arrays of holes, 
is based on work done on the representation of inclusions in a solid (three-dimensional) 
body. The displacement field in the exterior of such an inclusion, of which a cavity is a 
special case, can be expressed (Ref. 13) either in terms of the transformation strain in 
the cavity, as: 

= 7) b ' - 3 ‘ 31) 

or as: 

* = 16*7^1 - v ) Iv ^ 3 ' 32) 

in terms of the transformation stress, where 

fijk = (1 - 2i/)(6ijlk + Siklj) — SjkU + 3 Uljlk 

9ijk = (1 - 2 t/)(Sijlk + Siklj — Sjkli) 4* 3 IJjlk 

and r is the length of the vector / joining a volume element dv to the point of observation. 

The resulting stress field can be derived from the displacement field by differentiation. 
Evaluation of the displacement and stress fields by numerical integration is relatively 
straightforward. Routines for this evaluation can easily be designed to control accuracy 
based on the distance between the point at which the solution is to be calculated and the 
cavity. In particular, at locations remote from the cavity, the expression can be simplified 
to: 


Ux 


V _^J± _ 

16tt/x(1 — v)r~ 

v _^ 

8 tt (1 — i/)r 2 


In applying such solutions within the boundary element method, the transformation 
strain is chosen so that the sum of the transformation stress and the stress due to the 
boundary loads vanishes within the cavity. 
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Figure 3.2: Cube Containing a Cavity 


The decision to pursue this approach was guided by experience obtained during Task IIC, 
using a particular solution of this type. An exact (closed form) solution was obtained for 
an arbitrary, uniform transformation strain within a rectangular parallelepiped. Consid- 
erable experimentation was done to verify the practicality of the method outlined above. 
The test problem studied was that of a cubical cavity centered within a larger cube 
(Figure 3.2). The cube was subjected to an applied displacement along one axis. The 
transformation strain in the cuboid was defined by setting the integral of the total stress 
in the cuboid (using a 2x2x2 Gauss rule) to zero. The results obtained using the uniform 
transformation strain solution agreed well with results in which the cavity was modelled, 
over a wide range of cavity sizes. The load produced by the imposed deflection is shown 
in Figure 3.3. Even when the cavity occupied 50% of the total volume, the result of the 
approximate procedure produced a load only 8% lower than a complete modelling of the 
cavity. It is clear that the use of transformation strain solutions can correctly account for 
the overall load and displacement effects of a cavity, even at locations relatively near the 
cavity. Further, the results at relatively remote locations are quite insensitive to cavity 


26 



CAVITY SIZE, B/L 

— MODELLED + APPROX QmOD BEM ACUBE 


Figure 3.3: Load as a Function of Cavity Size 


shape. Figure 3.4 shows the comparison between the loads due to a cubical cavity and a 
spherical cavity of equal volume. The load is seen to be insensitive to cavity shape even 
at a cavity size ratio of .4, for which the load is computed only 1.5 cavity dimensions 
away from the cavity boundary. 

Further, even though the cuboid solution assumes uniform transformation stress, good re- 
sults are obtained for local variation of field quantities near the hole. Figure 3.5 compares 
results obtained for the x-displacement (Poisson contraction for the solid cube) both by 
directly modelling of the cavity and by use of the cuboid solution. In this case the cav- 
ity size ratio is .6 and the displacements are calculated on a line passing within .125 
cavity dimensions of the cavity. The presence of the cavity substantially modifies the 
Poisson contraction, which would be independent of axial location for a solid cube. The 
agreement between the results for the modelled and unmodelled cavities is excellent. 

The success demonstrated using the simple cuboid solution gives considerable confidence 
in the ability to model arrays of discontinuities using closed form or approximate inclusion 
solutions derived from (3.31) and (3.32). The technique will be computationally efficient 
and will fit within the same general framework as the particular solution formulation for 
inhomogeneity and thermal stress. 
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Figure 3.4: Effect of Cavity Shape on Load 

The only significant drawback of the method is the assumption of a uniform transforma- 
tion stress within each inclusion. This will become important only if the cavity is large 
enough that the stress varies significantly over the region occupied by the hole and if local 
results near the hole are required. In the short term this can be addressed by directly 
modelling a single hole in order to obtain the required local results. A more satisfactory 
long term solution would be the development of generalized versions of (3.31) and (3.32) 
allowing for variation of the transformation strain or stress. 

3. 3.2.7 GENERALIZED FUNCTIONAL REPRESENTATIONS The 
present version of the BEST3D program models part geometry using quadratic isopara- 
metric surface patches. Displacement and traction over each patch are represented using 
either quadratic or linear isoparametric shape functions. Although such representations 
can give adequate definition of both geometry and field quantities for almost any real 
problem, there are several undesirable features: 

1. The elements are not able to reproduce exactly a number of frequently used surface 
types, such as cylinders, spheres and ellipsoids; 
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Figure 3.5: Accuracy of Response Near Cavity 

2. The modeling of real parts using such elements requires considerable hand calcu- 
lation or the writing of special purpose computer programs, in order to fit surface 
patches to the true geometry; 

3. The elements are not entities used in the many geometric modeling systems which 
are, increasingly, being used for design purposes in the aerospace industries. 

A number of solid modeling systems presently under development, or already commer- 
cially available, use rational b-splines (in two dimensions) and rational b-spline surfaces 
(in three dimensions) for the representation of surface geometry. In such a representation 
each coordinate is represented in the form: 

S(s,t) = ( x(s 1 t),y(s,t),z(s,t )) (3.33) 

s)BUi)hii 

where the Pij and the hij are the control point values defining the geometry. The B h k are 
a suitable set of basis functions, normally taken to be piecewise cubic. These represen- 
tations possess the smoothness characteristic of real part geometry. Additionally, many 
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Figure 3.6: Geometric Model of Low Cycle Fatigue Specimen 


types of geometric details can be exactly represented using b-splines. It is, therefore, 
possible to provide a very accurate geometry representation using only a few elements. 
As an example, the low cycle fatigue specimen shown in Figure 3.6 consists of only 26 
rational b-spline surfaces. It was constructed in less than three hours using the CAEDS 
modeling system. Equally accurate representation of the same geometry using the current 
BEST3D modeling technology would require a minimum of 150 surface patches. 

The evident advantages of the b-spline representation have made it important to de- 
termine whether such a representation can be successfully incorporated in a boundary 
element analysis. In order to carry out this experimentation as rapidly and efficiently as 
possible, a two-dimensional (plane stress) boudary element stress analysis code using ra- 
tional b-splines to represent geometry, displacement and traction has been written. Key 
questions to be investigated include: 

1. The numerical stability of the resulting equation system; 

2. Sensitivity of the solution to the choice of collocation points; 

3. Whether or not displacement compatibility must be enforced as an additional con- 
straint at the break points between boundary segments. 



Figure 3.7: Beam Deflection (u = 0) 


The results to date are promising, although considerable work remains to be done. The 
test problem currently being used is that of a cantilever beam subject to a (parabolically 
distributed) shear load on the free end. This problem has always been an exceptionally 
difficult one for boundary element codes using isoparametric elements, but should have 
an essentially exact solution within the function space of rational b-splines. As is shown 
in Figure 3.7, when Poisson’s ratio is equal to zero, the spline based boundary element 
code gives excellent results over a wide range of beam aspect ratios. However, when 
Poisson’s ratio is nonzero (Figure 3.8), the results gradually deteriorate as beam aspect 
ratio increases. 

Investigation of the problem shows that the system equations are very poorly conditioned 
relative to the isoparametric formulation. The solution is extremely sensitive to the 
values of the off-diagonal terms in the traction kernel. Since the degree of sensitivity is 
the same for v — 0 and v — .35, the preliminary conclusion is that further attention is 
required to the accurate calculation of the traction kernel integrals. In addition, possible 
restructuring of the equation system to improve stability is being investigated. 
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Figure 3.8: Beam Deflection ( i / = .35) 

3.3.3 TRANSIENT AND DYNAMIC STRESS ANALYSIS 

3. 3. 3.1 REVIEW Three-dimensional problems of elastodynamics were not at- 
tempted until recently principally because of enormous computing requirements and the 
formidable task of numerical implementation. In order to reduce the computation and 
complications involved, simplifications of the BEM formulation dictated by the nature of 
the problem to be solved have been developed by a number of researchers. 

Dominguez (Ref. 17) simplified the steady-state dynamic kernel functions for the spe- 
cial case of periodic surface loading on rectangular foundations. He introduced simpler 
kernel functions to study the response of embedded rectangular foundations subjected to 
travelling waves. Karabalis and Beskos (Ref. 18) have done similar simplifications to the 
time domain transient boundary integral formulation. Yoshida, et. ai (Ref. 19) used a 
simplified BEM formulation for determining the response of a square foundation on an 
elastic halfspace, subjected to periodic loading and harmonic waves. Tanaka and Maeda 
(Ref. 20) have developed a Green’s function for two-layered visco-elastic medium, and 
using this Green function in a simplified BEM formulation they numerically calculated 
the compliances for a hemispherical foundation. More complex problems involving the 
periodic response of piles and pile groups have been attempted by Sen, et. al. (Refs. 
21-23) and Kaynia and Kausel (Ref. 24). They simplified the boundary integral formu- 
lation so that only displacement kernels are involved in the formulation. Some authors 
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have also introduced a potentially unstable method involving an ‘auxiliary boundary’ so 
that singular integration can be avoided. In all of the above works, the displacements 
and tractions are assumed to be constant over each surface element. 

Recently, Rizzo, eL al. (Ref. 25) and Kitahara and Nakagawa (Ref. 26) have imple- 
mented the BEM formulation for steady-state elastodynamic problems in a general form. 
Rizzo also implemented both a mixed-transform inversion to obtain the response in the 
time domain and a technique for the problem of fictitious eigenfrequency in certain ex- 
terior problems with homogeneous boundary conditions. Kitahara and Nakagawa have 
introduced series expansion of the periodic kernels for the low frequency range, to obtain 
a stable solution at low frequencies. 

In the present work, the direct boundary element formulations for periodic dynamic anal- 
ysis and time-domain transient analysis have been implemented for problems involving 
isotropic, piecewise homogeneous (substructured), three-dimensional solids. These imple- 
mentations are general and complete in all respects. In addition, for nonlinear, transient, 
dynamic analysis of three-dimensional solids, the direct boundary element formulation 
and its numerical implementation are presented for the first time. 


3. 3. 3.2 FURTHER DEVELOPMENT OF THE ELASTODYNAMIC FOR- 
MULATION During the course of Task IVC several new features have been added to 
the dynamic stress analysis part of the code to improve results as well as the efficiency 
of numerical computations. 


3. 3. 3. 2.1 Built-in Symmetry Capabilities In obtaining numerical solutions, a 
built-in symmetry capability allows one to solve problems having geometric and load- 
ing symmetry by modeling only a part of the actual geometry. The major steps in this 
procedure are briefly explained as follows: If the geometry and the boundary condition 
are symmetric with respect to a plane (or a number of planes), then only that portion of 
the boundary which lies on one side of the plane (or planes) is modeled. The symmetry 
can be with respect to y — z plane (half-symmetry), y — z and x — z planes (quadrantal 
symmetry), or y — z, x — z and x — y planes (octant symmetry). The effect of the unmod- 
eled part of the boundary is included according to the following scheme: For all the field 
points, the contributions of the unmodeled portions to the matrices of coefficients Fij and 
Gij are accounted for by reflecting the modeled surface elements with respect to the plane 
(or planes) of symmetry and then integrating over the reflected elements (with proper 
normals). For source points on the plane (or planes) of symmetry the contributions are 
added up directly whereas for all other source points the correct signs of the contributions 
are determined by the directions associated with the field variables with respect to the 
plane (or planes) of symmetry. By avoiding the calculation of identical quantities, this 
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procedure shortens the time required to evaluate the matrices. In addition, it reduces 
the time required to solve the set of linear equations, because the system matrix is much 
smaller. 


3. 3. 3. 2. 2 Linear Time Stepping and Interior Stress Calculations The direct 
boundary integral formulation for a general, transient, elastodynamic problem can be 
constructed by combining the fundamental point force solution of the governing equations 
(Stoke’s solution) with Graffi’s dynamic reciprocal theorem. Details of this construction 
can be found in Banerjee and Butterfield (Ref. 14). For zero initial conditions and zero 
body forces, the boundary integral formulation for transient elastodynamics reduces to: 

c tJ (Ou,ti,T) = [ [G t} (x,{,T) * t,(x,T) - F, } (x,tT) * u,(x,T )] dS(x) (3.34) 

where 

(3.35) 


(3.36) 

are Reimann convolution integrals and £ and z are the space positions of the receiver 
(field point) and the source (source point). The fundamental solutions G t j and Fij are 
the displacements and tractions at a point x and at a time T due to a unit force vector 
acting at a point £ at a time r. Equation (3.34) represents an exact formulation involving 
integration over the surface as well as the time history. It should also be noted that this 
is an implicit time-domain formulation because the response at time T is calculated by 
taking into account the history of surface tractions and displacements up to and includ- 
ing the time T. Furthermore, equation (3.34) is valid for both regular and unbounded 
domains. 

Once the boundary solution is obtained, the stresses at the boundary nodes can be 
calculated without any integration by using the scheme described for the static case. For 
calculating displacements at interior points equation (3.34) can be used with Cij = <5 tJ 
and the interior stresses can be obtained from 

= J[G' jk (x,£,T) * U(x,T) - F' k (x,£,T) * Ui(x,T)\ dS(x) (3.37) 

The functions GT and FT in the above equation are too complex to be listed here. 

In Task IVC constant time stepping was used to obtain the transient dynamic response. It 
was found to be more efficient to use a linear time variation of u and t on the boundaries. 


r T 

Gij*U = / Gij(x,T;£,r) U(x,t) dr 

Jo 

r T 

Fij *ui = / F tJ (x, T ; £, r) tq(z, r) dr 

Jo ~ 
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In this case: 


(3.38) 


N _ _ 

Ui(s.,r) = ^2[M i + M 2 «"(x)j 


n~ 1 

N 


U{x,t) = ^[Mi f"" l (i,r) + M 2 i*{x)] 


n=l 


where Mi and M 2 are the time functions, and are of the form: 


Ml = 

M 2 = L i^n(T) 


(3.39) 


(3.40) 

(3.41) 


For illustration purposes, consider the boundary integral equation for the first time step, 


i.e. 


djUi((,Ti) - JjPxjU — FijUi] dS dr - 0 


(3.42) 


The time integration in equation 3.42 by utilizing 3.39 is done analytically. After the 
usual numerical integration and assembly process, the resulting system equation is of the 
form: 

[4][X‘] - $][r,i + Mt-n - [B, , ][r“] = 0 (3.43) 

where: 

• A and B are matrices related to the unknown and known field quantities, respec- 
tively: 

• X and Y are the vectors of unknown and known field quantities, respectively: 

• for X and Y the superscript denotes the time: 

• for A and B the superscript denotes the time step at which they are calculated, 
and the subscript denotes the local time nodes (1 or 2) during that time-stepping 
interval. 


Since all the unknowns at time T = 0 are assumed to be zero, equation 3.43 reduces to: 

W = [BOX 1 ] + [B[][Y 0 } (3.44) 


For second time step, the assembled system equation has the form 



As in the constant time variation scheme, only the matrices on the right hand side of 
equation 3.45 need be evaluated. However, one needs to integrate and assemble four 
matrices at each time step as compared to two in the case of constant time variation. 
This can be done with only a small increase in computing time by integrating all the 
kernels together and then assembling all the matrices together. Equation 3.45 can be 
rearranged such that: 

[ 4P 2 ] = [B'JY 2 ] - [A\ + 4][jr'] - [b; + B-][y'] + [Bjin (3.46) 

In the above equation, all the quantities on the right hand side are known. Therefore, 
the unknown vector X 2 at time T 2 can be obtained by solving the above equation. 

Thus, for the present case, the boundary integral equation 3.46 can be written in dis- 
cretized form as: 


[a;][x w ] - [B ;i[y"] = - + a;- 1 ]!*"-"*'] - [B- + sr , ][y''-" +I ] + [sriim 

n=2 

(3.47) 


or 


K][**] = [B\W l N ) + [R N ] 


(3.48) 


It is of interest to note that, if time interpolation functions Mi and M 2 are replaced by 
Mi = M 2 = 0.5<^ n (r), the time stepping scheme for linear variation can be used for the 
case of constant variation with averaging between the local time nodes. 


3. 3. 3. 3 EXTENSION TO NONLINEAR TRANSIENT DYNAMICS In this 

section, a direct boundary element formulation and its numerical implementation for 
nonlinear, transient, dynamic analysis of three-dimensional, deformable solids of arbi- 
trary shape and connectivity is presented. The formulation is based on an initial stress 
approach, and is the first of its type in the field of boundary element analysis. The non- 
linearity considered in this analysis is that due to the nonlinear constitutive relations, i.e. 
material nonlinearity. The boundary integral equations are cast in an incremental form, 
and thus, elasto-plastic relations of the incremental type are used for material descrip- 
tion. These equations are solved by using a time-stepping algorithm in conjunction with 
an iterative solution scheme to satisfy the constitutive relations. The resulting algorithm 
is an unconditionally stable implicit scheme. However, the size of the time step that can 
be used is restricted by the size of the elements used for modelling the surface of the 
problem under consideration. 

The direct boundary integral formulation for a nonlinear, transient, dynamic problem, 
based on an initial stress approach, can be constructed by following a procedure similar 
to the one that has been used for a nonlinear static problem (Ref. 14). Under zero initial 


36 


conditions and zero body forces, the boundary integral equation for nonlinear transient 
dynamics is of the form: 

c^(e.T) = / s [G, i ( £ ,|,T)*<,(x,T)-F, J ( i ,e,T)*u 1 (x,T)]d5( £ ) 

+ / B ih (x, l, T ) * <r° ( (x, T) dV(x) (3.49) 

J v 

where * denotes convolution, and 

• f and x are the space positions of the receiver (field point) and the source (source 
point), respectively: 

• a°. is the initial stress tensor: 

• V denotes the interior of the body: 

• the fundamental solutions G^, Fij and Buj are displacements, tractions and strains 
due to an impulse in an infinite solid. 

Assuming all the field quantities to have a zero value at time T = 0, the boundary integral 
equation 3.49 can be written in an incremental form as follows: 

= f [G tJ (z, £,T) * A t t (x,T) - F tJ (x,^T) * A «,-(*, T)] dS(x) 

+ / 5 l7i (x,e,T)*A^,T)dF(x) (3.50) 

j v 

where A denotes the incremental quantity. 

The stress increment at an interior point 

* A«i(aE,T)] 

+ J v * A^foT) dV + Jihk^° u (i,T) (3.51) 

can be obtained by taking derivatives of equation 3.49 and using the constitutive rela- 
tionships 

Acr;y = DijkiAtid ~ A(J®. 

The functions G^ k , F* k , B*- k and Jujk are similar to their static counterparts except that 
some of the terms involve ( t — r/cp) and (t — r/c,), the retarded times for pressure and 
shear waves, respectively. 
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In (3.51), the volume integral must be evaluated in the sense of ( V — V € ) as V € — ► 0 as 
discussed in section 3. 3. 2. 2. The tensor Jujk is the jump term derived from the analytical 
treatment of the integral over V € . This jump term is the same as that of static plasticity 
and is independent of the size of the exclusion V €1 provided the initial stress distribution 
is locally homogeneous. 

The equations for incremental stresses cannot be constructed at the boundary points 
by taking the field point, f, in equation 3.51 to the surface, due to the strongly singu- 
lar nature of the integrals involved. However, the equations for incremental stresses at 
boundary points can be constructed by using a scheme similar to that described for the 
quasi-static nonlinear problem. Thus the stresses at a boundary point can be obtained 
froixu __ _ 

A *ik(£,T) = G' jk At,(i b ,T) -F^A^T) - B'. k A<T°(£ 6 ,T) (3.52) 

It should be noted that the above equation is free of any integration and time convolution. 

In dynamic plasticity, the choice of an appropriate constitutive model depends largely 
on the material properties and the loading conditions of the problem in hand. For this 
reason various constitutive models have been used for dynamic plasticity. For simplicity 
in the present analysis, the Von Mises model with isotropic variable hardening is used. 
The use of any other constitutive equations within the present formulation will not cause 
any difficulty, as long as the evoluton of A<r° can be obtained. 
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3.4 


COMPUTER PROGRAM DEVELOPMENT 


3.4.1 INTRODUCTION 

The computer program BEST3D was developed during Tasks IC and IIC to provide a tool 
for applying the three-dimensional boundary element method to the structural analysis 
of gas turbine engine hot section components. The program, described in detail in the 
First (NASA CR-174700) and Second (NASA CR-175060) Annual Status Reports, was 
designed to accommodate structures with very general geometry and loading. Further it 
was clear that additional capabilities would be developed for BEST3D over a period of 
several years. For this reason the program was designed so that the anticipated capabili- 
ties could be incorporated within the original framework and data structure of BEST3D, 
without requiring major recoding of existing capabilities. 

During Task I VC a number of new capabilities have been incorporated in BEST3D. 
The basic structure of the program remains, however, very similar to that of the original 
version. The major changes and additions made to the code during Task IVC, or presently 
in progress, are described in the following sections. 


3.4.2 GLOBAL PROGRAM STRUCTURE 

The major changes in the overall structure of BEST3D are: 


• The incorporation, in the static analysis, of a branch for natural frequency and 
mode shape calculations. 

• Addition of an additional nonlinear solution sequence (variable stiffness plasticity). 

• The implementation, in a test version of BEST3D, of a nonlinear dynamics capa- 
bility. 


The overall structure of BEST3D at the end of Task IVC is shown in Figure 3.9. The 
program consists of a common input section, followed by three branches, for static, forced 
response and transient analysis. The static analysis branch is the model for the entire 
code, since the other branches largely employ generalized forms of the same algorithms 
used in the static analysis. The branch used for natural frequency /mode shape calculation 
is actually part of the static analysis loop. 
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Figure 3.9: BEST3D - Overall Program Structure 














3.4.3 PROGRAM INPUT 


The input to BEST3D is essentially unchanged. Input formats are described in the 
BEST3D User’s Manual (Ref. 15). The minor input changes made during Task IVC in- 
volve the addition of parameters to control the nonlinear dynamics and natural frequency 
calculations. 

The BEST3D User’s Manual (Ref. 15) is periodically updated to remain current with 
BEST3D development activity. The current version of the manual is consistent with the 
code presently available at the Lewis Research Center of the National Aeronautics and 
Space Administration. 

3.4.4 IMPLEMENTATION OF VARIABLE STIFFNESS PLASTICITY 

Equations 3.11 and 3.12 can be discretized and assembled into the following sets of 
equations by separating the known and unknown boundary components: 

A b x = B b y + C b K\ 

u = A i i + B?y + &K\ (3.53) 

A = LA'x + LB'y + LC"K\ 

where the matrices A, B, C are constant for a given problem and the matrices AT, L which 
are dependent upon the state variable, are assumed to be constant during each small 
load step. The superscripts b and i denote boundary and interior related quantities, 
respectively. 

It should be noted that the first equation is written for the boundary nodes only, the 
second only for the interior nodal points (excluding boundary nodes), and the third 
equation is written for all (boundary as well as interior) nodes likely to yield as a result 
of any applied loading. 

The first and last of the above equations can be rewritten as: 



A b = b b + C Xb A 

(3.54) 

and 

A = A x i + b x + C x \ 

(3.55) 

or 

II 

>- 

+ 

O-* 

>• 

(3.56) 

where 




b b = B b y , C Xb = C b K , A x = LA ' 
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b x = LB'y, C x = LC'K, H = I - C x 


and I is the identity matrix. 

Equation 3.55 can be recast as: 

A = A°x + 6° 

where 

A 0 = H~ X A X 
b° = H~ l 2 3 4 5 6 b x 

Substituting the above equation into 3.54 results in the final system equations: 

A r i = b r (3.57) 

where 

A r = A b - C Xb A° 

b r = b i -C x 

The above equation can be solved to evaluate the unknown vector x at boundary nodes 
for every increment of loading. It should be noted that the matrix A b is the reduced 
system matrix and b b the corresponding right hand side vector of a similarly posed elastic 
problem. The present formulation is similar to the variable stiffness approach used in the 
finite element method since the system matrix on the boundary as well as the right hand 
side vector is modified for each increment of loading. 

The solution process does not involve any iterative procedure, instead the substantial 
part of the solution effort is spent on assembly of the system equations for each load 
step. These operations can be described as follows: 

1. Impose an arbitrary boundary loading and solve the elastic problem in the usual 
manner. 

2. Scale the elastic solution such that the highest stressed nodes are at yield. 

3. Apply a small load increment and compute K and L matrices using the past history. 

4. Form the system equation and solve for x from 

A r i = b r 


5. Evaluate A from 


A = A° x + b° 


6. Evaluate & and u and accumulate the incremental quantities. 
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7. Return to step (3) if the displacements are less than a specified norm. Otherwise, 
failure is assumed to have occurred. 

It is of interest to note that the matrices if, I do not exist in the elastic region, therefore, 
the corresponding equations involving these matrices are formed, and A is determined, 
only for the plastic nodes. 

Further, the load step must be small (usually less than five percent of yield load). When 
the yield condition is violated, an initial stress equilibrium correction is applied to bring 
back the stress point to the yield surface and this correction is taken forward to the next 
load step. 


3.4.5 INCORPORATION OF PARTICULAR INTEGRALS 

The branch designed for the incorporation of particular integrals for the solution of prob- 
lems with thermal loading, inhomogeneity and/or embedded holes or cracks is shown in 
Figure 3.10. It is closely related to the algorithm, described in the Second Annual Status 
Report (NASA CR-175G60), for the determination of natural frequencies. Following the 
calculation and storage of the coefficient matrices for both the boundary and interior 
stress equations, the displacements and tractions due to the particular solution are calcu- 
lated. The full matrices containing the particular solution values are never stored, since 
the required multiplications with already existing boundary element matrices are carried 
out as the calculation proceeds. After the system matrix assembly the weights used in 
the initial strain approximation are eliminated from the system, leading to a modified 
system matrix similar to that used in variable stiffness plasticity. This matrix is then 
decomposed and the remainder of the problem solved exactly as in the original code. 


3.4.6 EXTRACTION OF EIGENVALUES 

The routine used for the solution of the generalized algebraic eigenvalue problem is based 
on the work of Ref. 16. The original paper addressed the extraction of the largest few 
eigenvalues of a very large, very sparse system arising in the development of multi-grid 
methods. Two processes, an iteration and a purification step, were used in the original 
method. To date, only the iterative process has been required in the BEST3D application. 
In order to employ the algorithm in BEST3D it was necessary to reformulate it for the 
generalized eigenvalue problem and adapt it to the block storage used in BEST3D. 

The generalized eigenvalue problem: 


(A — u 2 B)x ~ 0 
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Figure 3.10: Incorporation of Particular Integrals 


in which the lowest eigenvalues are sought, is reformulated as the ordinary eigenvalue 
problem: 

(A~ l B- \I)x = 0 

with, 

A — 1 / aT 

in which the largest few eigenvalues are sought. Starting from an arbitrary initial vector, 
x°, successive vectors orthogonal to one another with respect to the matrix A 1 B are 
formed. At each step of the iteration an upper Hessenberg matrix is formed, whose 
eigenvalues converge to the eigenvalues of the original system. Convergence is excellent, 
with the first mode usually obtained in five to eight iterations. Computational experience 
has shown that one additional mode per iteration is normally obtained thereafter. 
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3.4.7 NONLINEAR DYNAMICS 


Equations (3.50) and (3.51) provide the formal basis for developing the dynamic plasticity 
algorithm. However, the initial stresses Acr^. defined in equations (3.50) and (3.51) are 
not known a priori and have to be determined by satisfying the constitutive relations. 
Thus, equations (3.50) and (3.51) and (3.52) can be regarded as a coupled system of 
nonlinear equations. In the present implementation, equation (3.51) and (3.52) are used 
to calculate the stresses at interior and boundary points and the nonlinear material model 
is then used to evaluate the inelastic stresses. Since the volume integrals of inelastic stress 
vanish except in regions of nonlinear material response, approximations of geometry and 
field quantities are required only where nonlinearity is expected. In the present work, 
isoparametric (quadratic) volume cells are used for approximating the geometry and the 
variation of initial stresses such that: 

Xi = Mp^rfjzip (3.58) 

< = 

where z, are cartesian coordinates, 2,73 are nodal coordinates of the volume cell, Mp (f3 = 
1,20) are the quadratic shape functions for the volume cell, and the overbar denotes nodal 
quantities. 

The volume integral of equation (3.50) can then be represented as: 

J q T J v B tlk ( X, T; i b , r) Ac^z, r) dV dr = (3.59) 

£ [ T / £<y(* m (2),T;| 4 ,T) M 0 (r,) dV m A* 0 ” dr 

m=l J0 JV ’ n 

where: 

• £ is the field point on the boundary (boundary node), 

• 2L m (v) is a point in cell m, 

• Aov™ are the nodal values of incremented initial stress for the m th cell, 

• V m is the m ih volume cell, and 

• L is the total number of cells in a single region. 

Similarly, the volume integral of the interior stress equation (3.51 ) can be expressed as: 

/. r Jr S .y T; £’ T > T ) iV iT = ( 3 ' 6 °) 
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£ f* / B' ,(*», T;£, r) M„(j) <iV m A<r£ <ir 
in which the time integral is treated analytically as before. 

The volume integration scheme used in this work is substantially similar to the quasi- 
static nonlinear analysis except that the cells are sub-divided to account for the silent 
part of the cell which has not yet been affected by the shear and pressure waves. 

3.4.7. 1 TIME-STEPPING In order to obtain the nonlinear transient response at 
a time T„, the time axis is discretized into N equal time intervals, i.e. 

N 

T n = £ nAT 

n=l 

where AT is the time step. 

Using equation (3. 4. 7.1), the integral equation (3.50) can be written as: 

c^At*;(f,Tjv) — / f (GijAti — FijA-Ui) dS dr 
JTh.x Js 

- T* ( B t u Ac., dV dr = (3.61) 

JTs-i JV 

[ " ' [ (G^Ati - FijAiLi)dS dr + f BujAcr^dVdr 
J r=0 Js Jr — 0 

For the present case the linear time interpolation scheme is used to approximate the time 
variation of the field quantities during a time step because the same scheme can also be 
used for constant time interpolation with averaging. 

After the usual discretization and integrations (in time and space), the integral equations 
(3.61) are transformed into an assembled system equation of the form: 

[^[AX"] - [Bj][Ar w ) - [c;][a/'"i 



= - E [i A 2 + A n 1 - 1 }[AX N ~ n+1 } - [J r 2 + Bp l ][AF w - n+1 ] 

n=2 

(3.62) 


+ [C 2 n + C 1 n - 1 ][Aa°’ JV - n+l ]] 


or 

[a\][ax n ] = {b\](ay n ] + [c!][aO + [a?] 

(3.63) 

or 

A b AX N = A b bN + A b° M 

(3.64) 

where: 
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• A and B are the matrices related to the unknown and known incremental displace- 
ments and tractions; 

• C is the matrix related to the initial stresses; 

• AX and AY are the vectors of unknown and known incremental displacements and 
tractions; 

• for AX, AY and Act 0 , superscripts denote time, i.e., AX n = X n — I”" 1 ; 

• for the A, B and C matrices, superscripts denote the time step when they are 
calculated, and subscripts denote the local time node (1 or 2); 

• R N is the effect of past dynamic history; 

• and, 

A b bN = [B*][AY*] + [R N ] 

A b - 

A b°' bN = [Z?J][Aa 0>JV ] 

Similarly, the integral equation for stresses can be written in a discretized form as: 

[Aa"] = [ijKAjr"] + [Sj][Ar"] + [c^a/'"] + [jr"] (J.«5) 

or 

Aa N = A"AX N + A b* N + A b°'* N (3.66) 

where the overbar indicates that the matrices are related to the stress equation: 

a* = tQ; = [bI][ay n ] + [r* n ] 

and; 

A b 0 '*" = [C^][Acr 0 ’"] 

3.4. 7. 2 ITERATIVE SOLUTION ALGORITHM FOR DYNAMIC PLAS- 
TICIT Y The algorithm described here provides the solution of system equations (3.62) 
and (3.65). The solution of these system equations requires complete knowledge of the 
initial stress distribution within the yielded region that is induced by the imposition of 
the current increment of boundary loading. This, unfortunately, is not known a priori for 
a particular load increment and, therefore, an iteratve process must be employed within 
each time step. 

This incremental algorithm can be described as follows: 
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1. Obtain the transient elastic solution for an arbitrary increment of boundary loading 
A Y N during the time interval T„_i to T„, a s 

A b AX N = Ab bN 

and 

a/ = A* AX N + Ab' N 

where N is the time step number. If the material has not yet yielded, accumulate 
X- vectors, i.e., X N = X* -1 + AX N . 

2. If the material has previously yielded, go to step (6). 

3. Check whether any node has yielded during the current time step. If the material 
has not yielded yet, accumulate stress and strain, and go back to step (1). 

4. Calculate the value of cr°, the equivalent stress, by using cr T = cr N ~ l + Acr N as the 

stress changes and compile a list of yielded nodes. For elastic nodes accumulate the 
stress and strain, i.e., cr n = cr T and + [D*)~ l Acr N . Calculate the correct 

stress at the elastoplastic nodes by using the elastoplastic stress-strain relations 
A cr** = D ** Ae and using the elastic strain increments as a first approximation. 
Modify the stress history for yielded cells cr N = aN — 1 + Acr^. Calculate initial 
stress, A a — a — a . 

5. Assume A6^ = 0 and A b bN — 0 and, using the generated initial stress Act 0 , 
calculate a new AX N by using equation (3.64) and A <j N by using equation (3.66). 
Calculate the equivalent stresses by using the history a T = cr^-h Acr M and compile a 
list of yielded nodes. For elastic nodes, accumulate the stress, <j N = cr T , and strain. 
For the elasto-plastic nodes calculate the current stress increment, A<r* p . The initial 
stresses generated are Acr° = Acr N — A cr* 1 *. Modify the stress history for the yielded 
nodes, cr N = cr N + A^. Accumulate AX N and A<j°’ N (i.e.= AX^ = &X N + AX N 
and A(7 0iV = Acr 0 ^ + Act 0 ) so that they can be used in the next time step for past 
convolution. 

6. Check if the initial stresses Act 0 are less than an acceptable norm. If so, go to step 
(1). If not, go back to step (5). If the number of iterations exceeds a fixed number 
(usually 50), it is assumed that collapse has occurred. 


48 



3.5 


VALIDATION/VERIFICATION 


3.5.1 SUMMARY 

During Task IVC, significant new capabilities were added to BEST3D. In addition, new 
algorithms presently under development were tested using two-dimensional programs or 
special purpose three-dimensional codes. The validation and verification of these new 
capabilities and the continuing general verification of the BEST3D code are discussed 
in the subsections which follow. Attention in this report is directed primarily to the 
dynamic capabilities of BEST3D, the initial evaluation of variable stiffness plasticity and 
the new particular integral formulations and the continuing application of BEST3D to 
three-dimensional component analysis. 

Table 3.1: Calibration of Major BEST3D Capabilities 


BEST3D CAPABILITY 

RESULTS USED 


EXACT 

OTHER NUMERICAL 

EXPERIMENTAL 

ELASTICITY 




general capabilities 

X 



stress concentrations 

X 

X 

Benchmark Notch 




Diffuser Case Boss 

thermal stress 

X 



PLASTICITY 




general capabilities 

X 

Finite Element 




Alternate BEM Algorithms 


structural details 

X 

Finite Element Analyses 

Theocaris /Marketos 



for Experimental Programs 

Benchmark Notch 

thermoplasticity 

X 


Carmen/ Hess 

DYNAMICS 




time domain 

X 

Pao/Mow 


frequency domain 

X 


Woods 

EIGENVALUES 

X 

X 

Keilb, et. ai 


The major capabilities of BEST3D together with analytical and/or experimental used to 
calibrate the capabilities are summarized in Table 3.1. 
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3.5.2 MULTI-REGION INTERFACES 


One of the key features of the BEST3D code is the ability to join two subregions across 
one or more interface elements, with either complete fixity or tangential sliding allowed 
at the interface. This capability is a key to the successful modelling of complex geometry 
and of the interaction between parts such as blades and disks. In order to verify this 
capability in BEST3D, a highly idealized model of one-half of a compressor dovetail was 
built. The model, shown in Figure 3.11, consists of twenty eight quadratic elements. The 
model was held on rollers and loaded with a radial load on the simulated blade, with 
centrifugal load and with combined centrifugal and radial loads. Cases were run with 
both fixed and sliding interface conditions. 

In order to check the correct operation of the BEST3D code for these problems, for 
which no analytical solutions are available, cases were run in which the modulus of the 
simulated disk (Region 1) was set to be several orders of magnitude larger than that in 
the simulated attachment (Region 2). In such cases the response in the attachment is 
the same as that obtained when the attachment alone is run with either fixed or sliding 
contact on the interface. 

For all of the cases run, the results of the single region analysis agreed with the two 
region, multi-material analysis. 


3.5.3 PLASTICITY USING VARIABLE STIFFNESS ALGORITHM 

The accuracy and efficiency of the newly incorporated variable stiffness plasticity algo- 
rithm in BEST3D has been investigated using the problem of the radial expansion of a 
thick, internally pressurized cylinder. The ratio of outer to inner radius of the cylinder is 
2.0. The material has a Young’s modulus of 2600, Poisson’s ratio of 0.3 and a yield stress 
of 600. The piecewise linear approximation to the stress-strain curve is shown in Figure 
3.12. The BEST3D model consisted of eighteen boundary elements and four (twenty 
node) volume cells. 

Three analyses were run: 

1. iterative algorithm with load increments of 5 percent of the yield load for the first 
seven increments, with 2 percent load increments thereafter; 

2. variable stiffness algorithm with 2 percent load increments throughout; 

3. variable stiffness algorithm with load increments identical to those used for the 
iterative algorithm. 
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EQUIVALENT STRESS 



Figure 3.11: Simulated Disk/Attachment BEST3D Model 



EQUIVALENT STRAIN 
Figure 3.12: Stress-Strain Response 


5 





Figure 3.13: Load - Displacement Response 


All three analyses gave virtually identical load-displacement curves, at both the inner 
and outer radii (Figure 3.13). Further, there is excellent agreement in the stress variation 
through the thickness of the cylinder (Figure 3.14), and in the tracking of the stress-strain 
curve (Figure 3.12). 

The run time for the third (variable stiffness) analysis (on the SUNY-B HP9000) was 1.6 
times that required for the iterative analysis (analysis 1). The current results indicate 
the possibility, presently being investigated, of using even larger load increments for the 
variable stiffness analysis. This would lead to considerable reductions in computing time, 
since variable stiffness computing time is directly proportional to the number of load 
increments. Further, very significant opportunities exist for improving the efficiency of 
the variable stiffness algorithm, as compared to the iterative algorithm, which is already 
the product of ten years of development. 


3.5.4 THERMAL/INELASTIC ANALYSIS USING PARTICULAR INTE- 
GRALS 

The major limitation in use of the boundary element method as a general purpose analysis 
tool is that a volume discretization and integration must be introduced to account for 
inelastic response, elastic inhomogeneity and/or non-steady state thermal strains. This 
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Figure 3.14: Stress Variation through Thickness of Cylinder 


requirement for volume discretization meant that the boundary element method became 
less attractive relative to the finite element method, in which inhomogeneity and thermal 
strains can be incorporated within the same framework as homogeneous elastic analysis. 

A high priority in the current phase of the inelastic methods program has been the 
reduction (or, preferably, elimination) of volume integration requirements. An improved 
analytical formulation has been gradually developed, starting from the exact initial strain 
solution for a cuboid and the development of particular integrals for centrifugal loads. 
The development of a formulation allowing solution of thermal stress problems and of 
problems involving smoothly varying elastic inhomogeneity (such as that caused by the 
variation of Young’s modulus with temperature) has now been completed. 

In order to evaluate the potential benefits of this new approach the simplest case, the 
particular integral for thermal strains, has been incorporated in BEST3D. Since thermal 
stress capability using volume integration is already provided in BEST3D, direct compar- 
ison is possible between the new and old algorithms. The test problem shown in Figure 
3.15, is geometrically simple, but very difficult to solve accurately. The difficulty arises 
from the fact that all the loading in the problem is generated by volume integration, 
since there are no mechanical loads applied. As a result the solution is very sensitive to 
volume discretization (especially cell aspect ratios) and to volume integration accuracy. 
In this analysis the cylinder is allowed to grow on rollers in response to a uniform tem- 
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ANALYTIC SOLUTION, PLANE STRESS CONDITIONS 
RADIAL DISPLACEMENT u r = a (AT)r, r = RADIAL DISTANCE 
STRESS a jj = 0.0 

Figure 3.15: Boundary Element Model of a Thin Disk 

perature increase. The analytical solution to this problem is a linear variation of radial 
displacement with radius and a zero stress state. Two cases were analyzed. In the first 
the disk thickness was .5 inch and in the second it was 0.1 inch. These cases produce cell 
aspect ratios of 5:1 and 25:1, respectively. The results are summarized in Tables 3.2 and 
3.3. 

It can be seen that the particular integral approach, for both cases, produces much more 
accurate results than the volume integral approach. Even in the most demanding case, 
the worst errors are .2% in displacement and 1.2% in stress (relative to Ea(AT)/(l — 2i/)). 
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Table 3.2: Free Thermal Expansion of a Thin Disk — (b — a)/t = 20 


■ 

Type of 
Analysis 

Nodes 

i 

2 

3 

4 

5 

6 


Volume Integration 

- 

- 

- 

9.41 

14.272 

18.872 

Ur 

Particular Integrals 

- 

- 

- 

9.99 

14.98 

19.980 


Analytical Solution 

- 

- 

- 

10.0 

15.0 

20.0 

<Jyy 

Volume Integration 

983.8 

1113.5 

1119.2 

218.9 

-177.8 

-346.9 


Particular Integrals 

1.6 

6.1 

3.1 

-5.9 

0.7 

-5.4 


Table 3.3: Free Thermal Expansion of a Thin Disk — (6 — a)/t = 100 



Type of 
Analysis 

Nodes 

1 

2 

3 

4 

5 

6 


Volume Integration 

- 

- 

- 

8.6 

11.5 

14.2 

Ur 

Particular Integrals 

- 

- 

- 

9.972 

14.964 

19.946 


Analytic Solution 

- 

- 

- 

10.0 

15.0 

20.0 

CTyy 

Volume Integration 

259.1 

405.5 

532.2 

-1104.9 

-1034.1 

-1494.6 


Particular Integrals 

-75.9 

-70.9 

-74.4 

-10.5 

-5.9 

-18.8 


Recently the problem of the inelastic expansion of a thick cylinder described in Section 
3.5.3 was analyzed by this newly developed method. The results were virtually indis- 
tinguishable from the iterative plasticity results already shown in Figures 3.12 through 


3.14. 










4-P(t) = 1000 


mm 

bbi 


■ 

a 


e 

E 


2 

s 

8 
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— ANALYTICAL 
• NUMERICAL (At = 0.00356) 
O NUMERICAL (At = 0.00445) 



2000.0 


1000.0 


0.0 


Figure 3.16: Dynamic Stress in Square Cantilever under End Load 


3.5.5 BAR SUBJECTED TO TRANSIENT END LOAD 
1. Square cross-section : 

A bar with square cross-section is held along its sides by lubricated rollers and is fixed 
at one end. The free end is subjected to a suddenly applied and maintained uniform 
compression t z = 1000. The dimensions of the bar are L = 8.0 and 6 = 2.0. In view of 
the material properties, the characteristic time required for the compressive wave to reach 
the fixed end is 0.0578 sec. Figure 3.16 shows the discretization and the numerical results 
for the normal stress, in which the results from the time domain algorithm for two differ- 
ent time steps AT are compared with the exact analytical solution for one-dimensional 
stress wave propagation. Although the numerical results are in good agreement with the 
analytical solution, it is clearly very difficult to reproduce the sharp jump in the stress 
as the disturbance reaches the point initially and when the reflected stress wave returns 
to the same location. This difficulty has been observed elsewhere as well. 

The axial displacement history at the free end is shown in Figure 3.17. The displacements 
are normalized by static displacements and the time is normalized with respect to the 
characteristic time required for the compressive wave to reach the fixed end. It can be 
seen that the numerical results are in good agreement with the analytical solution. The 
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1.50 



Figure 3.17: Deflection at Free End of Square Cantilever 


differences are mainly due to the three-dimensional nature of the simulated problem. 

2. Circular-<;ross-sectiQn : 

In order to investigate the effects of the cross-section on the numerical results, a bar 
with circular cross-section having the same material properties and boundary conditions 
as described in the last example was analyzed. The boundary element mesh for this 
problem is shown in Figure 3.18. The bar has a length L = 5 and diameter d = 1. Thus, 
the characteristic time required for the compressive wave to reach the fixed end is 0.02236 
sec. The time step used in this example is AT — 0.004475 sec. 

Figure 3.19 shows the numerical results for the normal stress at the midspan of the bar 
compared with a one-dimensional analytical solution. As mentioned in the last example, 
the sharp jumps in stress are diffused in the numerical results. However, by using more 
elements and smaller time steps, the numerical results in the vicinity of the jumps will 
agree more closely with the analytical solution. 

The time history of the normalized, axial displacements at the free end is plotted in 
Figure 3.20 against the one-dimensional analytical solution. The results are in good 
agreement, except for the peak displacements. The numerical peak values are less than 
than of the analytical solution and this results in an increase in the difference between 
the two solutions at later times. The difference, once again, is mainly due to the three- 
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Figure 3.18: Surface Discretization of a Circular Bar 


ANALYTICAL 

• BEM (AT = 0.004475s) 
p 0 = APPLIED PRESSURE 



TIME ( T x 0.004473 SEC) 


Figure 3.19: Axial Stress at Midspan of Circular Bar 



u z /u 


ANALYTICAL 

BEM (AT = 0.004475s) 





Figure 3.20: Axial Displacement at Free End of Circular Bar 
dimensional nature of the problem under consideration. 

3.5.6 SPHERICAL CAVITY SUBJECTED TO DYNAMIC PRESSURE 

A spherical cavity is embedded in an infinite medium with E — 8.993xl0~ 6 , v — 0.25, 
and p = 2.5xl0“ 1 * * 4 . The radius of the cavity is R = 212. Using the built-in symmetry 
capabilities, this problem is modeled by one octant only. The single octant modelled 
is idealized using three quadratic surface patches. The characteristic times required for 
the pressure and shear waves to travel a cavity radius are 0.00102 sec and 0.00277 sec, 
respectively. Four cases are considered: 

1. Spherical cavity under sudden radial expansion 

A radial pressure p = 1000 is suddenly applied and maintained at the cavity sur- 

face. Figure 3.21 shows the time variation of deviatoric stress at the cavity surface 

obtained by the time domain algorithm. Concurrently plotted is the result reported 
by Hopkins (Ref. 27) based on the work of Hunter (Ref. 28). In general, the nu- 
merical results are in good agreement with the analytical solution. The transient, 
time-domain solution remains stable and reaches the expected static solution at 


59 




Figure 3.21: Deviatoric Stress at the Cavity Surface 


larger times. It can be seen that the maximum deviatoric stress for the transient 
case is 1.77 times the applied pressure, whereas for the static case it is 1.54 times 
the applied pressure. 

2. Spherical cavity subjected to a triangular pulse of radial pressure 

A triangular pulse of radial pressure, as shown in Figure 3.22, is applied at the 
cavity surface. This example is solved by using linear time interpolation functions 
and two different time steps. The radial displacements at the cavity surface are 
plotted in Figure 3.22. The numerical results from both the time steps are almost 
identical. Thus, this example once again demonstrates the stability of the present 
algorithm. 

3. Spherical cavity subjected to a rectangular pulse of radial pressure 
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0.010 



Figure 3.22: Radial Expansion of Cavity under Triangular Pulse 


A rectangular pulse of radial pressure, as shown in Figure 3.23, is applied at the 
cavity surface. This example is also solved by using linear time interpolation func- 
tions and two different time increments. Figure 3.23 shows the time history of the 
radial displacement of the cavity. By comparing these results with those due to a 
triangular pulse (i.e. Figure 3.22), it can be seen that, in general, displacements at 
any time interval due to the rectangular pulse are twice that due to the triangular 
pulse. This is because the response depends upon the total impulse, and the to- 
tal impulse due to the rectangular pulse is double that due to the triangular one. 
Hence, the displacement amplitude response due to the rectangular pulse is also 
approximately double that of the triangular pulse. 

4. Spherical cavity engulfed by a pressure wave 

A propagating plane pressure wave whose front is perpendicular to the z-axis (Fig- 
ure 3.24) first impinges on the pole with coordinates (0,0,212). The resulting non- 
zero incident stresses are cr* = —1000, <f = <j 1 = (v/(l — v))a . This wave 

propagation type of problem has been solved by superposition (Miklowitz, Ref. 
29). A mesh using three quadrilaterals per octant (Figure 3.25) is employed here in 
conjunction with the time-domain approach. Figure 3.26 plots the hoop stresses 
<7*^ and cr* 0 normalized by the magnitude of the incident stress a\ z versus the 
non-dimensional time r* = RT/c . 

The plots are for three locations on the surface of the cavity: the two poles (< p = 0, 7 r) 
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RADIAL DISPLACEMENT 



Figure 3.23: Radial Expansion of Cavity under Rectangular Pulse 

Pressure Wove 



Figure 3.24: Spherical Cavity Engulfed by a Pressure Wave 
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X 


x 2 


X, 

MODE AMO CLEMENT CONFIGURATION 
WITH 24 RECTANGULAR CLEMENTS 

Figure 3.25: Mesh Used for the Analysis of Spherical Cavity 

and the equator ( <j> = x/2). Concurrently plotted are the results (Pao and Mow, 
Ref. 30; Norwood and Miklowitz, Ref. 31), obtained by analytical inversion of 
the Fourier transformed solution. Good agreement is observed between the two 
solutions. Finally, Figure 3.27 plots the radial displacement time history at the 
same three locations. 


3.5.7 DYNAMIC LOADING ON A SQUARE FOOTING ON A HALF 
SPACE 

In this example, a square flexible footing on a half-space (the soil under the footing) is 
subjected to time dependent vertical tractions. The mesh for this problem is shown in 
Figure 3.28(a). The side of the footing is B = 26 = 2, and the material properties of the 
half-space in any consistent units are: elastic modulus E = 2.6, Poisson’s ratio = 0.3 and 
mass density = 1.0. 
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Figure 3.26: Hoop Stress Induced by Passage of Pressure Wave 
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Figure 3.27: Radial Displacements Induced by Passage of Pressure Wave 
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Figure 3.28: BEST3D Models for Analysis of Footing Compliance 


The foundation is first subjected to a uniform harmonic vertical displacement u z of ampli- 
tude equal to unity. The surface of the halfspace is traction free. The traction distribution 
under the foundation obtained by the boundary element method is integrated to give the 
total vertical load P £ . The normalized compliance of the foundation in the vertical di- 
rection is obtained as C vv — fibu l /P z . T 

modelling the foundation as well as surface of the halfspace. Since the transtorm domain 
computer program can take advantage of symmetry, only one-quarter of the problem 
needs to be discretized. The coarse mesh uses four and twelve elements to model the 
foundation and the halfspace, respectively. Note that the outermost four elements are 
infinite elements. This discretization results in 44 nodes. The finer mesh uses six and 
twelve elements for the same purpose, leading to two infinite elements and 65 nodes. 

This problem wets originally solved by Wong and Luco (1976). They numerically inte- 
grated the vertical displacement at the surface of a homogeneous halfspace due to a unit 
point load over the foundation, which was discretized into small squares. This problem 
was recently revisited by Rizzo, eL ai (1985) using a BEM approach. In their work 
both frictionless and welded cases are considered and two approaches are used: The ex- 
act one employs the halfspace kernels (Lamb’s solution) and the approximate one uses 
the fullspace kernels (Stoke’s solution). In both cases, only the rigid foundation is dis- 
cretized and these two approaches are practically indistinguishable except at the very 
low frequency range. All three solutions mentioned are plotted in Fig. 3.29, along with 
the vertical compliance obtained by the present method, using the fine mesh. The good 
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Figure 3.29: Vertical Compliance of Flexible Square Footing 


agreement between the present results and that of Rizzo, et. al. (1985) should be no- 
ticed. The major difference between Wong and Luco’s results and the boundary element 
results is due to the fact that quadratic shape functions are used for representation of the 
variation in the field variables over each element in the present work, as well as that of 
Rizzo, whereas Wong and Luco assumed that the unknown contact stresses are uniform 
within each element. Differences in the results from the coarse and fine meshes are quite 
small (Table 3.4). 

Next the case of a transient loading was considered. The time step used for this analysis 
is AT = 0.2. The time history of the applied pressure and the vertical displacements at 
the center and corner of the loaded area are plotted in Figure 3.30. It can be seen that 
the vertical displacement at the center of the footing converges to the static value after 
2.4 seconds, whereas the vertical displacement at the corner of the footing converges to 
its static value at a later time. The mesh used for this problem gives a maximum error 
of 2% for static analysis. The results obtained for the present dynamic problem should 
have equivalent accuracy. This example shows the utility of the present algorithm for 
transient dynamic analysis of halfspace problems. 
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Table 3.4: Comparison of Vertical Compliances Obtained Using Two Different Meshes 



3?[CU(ao)] 


a 0 = ujb/c 2 

Coarse Mesh 



Fine Mesh 

0.5 

.118 

.117 


■BUS 

1.0 

.064 

.069 


mSSm 

1.5 

.032 

.034 


-.070 

2.0 

.021 

,018 

-.059 

-.052 

2.5 

M 

.015 

-.052 


3.0 

Km 

.012 

-.036 


3.5 

I 

.006 

-.035 

-.032 

4.0 

.004 

.004 

-.027 

-.027 


3.5.8 NONLINEAR DYNAMIC RESPONSE OF A BAR 

A bar with circular cross-section (of Section 3.5.5) is held along its sides by lubricated 
rollers and is fixed at one end. The free end is subjected to a suddenly applied and 
maintained uniform compression t £ = —333, which exceeds the yield stress of the bar 
(yield stress = 300). In this example, the bar has dimensions and material properties 
identical to that of the example discussed in Section 3.5.5. The discretization of the 
bar is similar to that shown in Figure 3.18 except that, in the present example, the 
full cross-section of the bar is modeled instead of one-quarter of it. The volume of the 
bar is discretized by using five equal twenty node volume cells. A bilinear stress-strain 
relation, as shown in Figure 3.31, is assumed. The time step used for this example is 
AT = 0.004473. Use of a uniform time step is vital, as it allows reuse at each new time 
step of previously calculated matrices. 

In Figure 3.31, the elasto-plastic response of the bar at timeT = 0. ST e (where T e = ciT/Jt, 
i.e. the time taken by the compression wave to reach the fixed end of the bar) is plot- 
ted against the one-dimensional analytical solution (Garnet and Armen, Ref. 32). The 
normal stress a LZ is normalized by the elastic modulus and the distance along the bar is 
normalized by the length of the bar. The numerical results are in reasonable agreement 
with the analytical solution except for the sharp jumps in the stress which are diffused by 
the numerical analysis. The analytical solution of course can never match the real mate- 
rial response where it behaves more like the numerical solution. The major differences in 
the results between the two solutions can be attributed to the three-dimensional nature 
of the present example. As the bar is on lubricated rollers, in addition to longitudinal 
stress, lateral stresses also exist in the bar. Simple one-dimensional theory considers 
longitudinal stress only. 
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Figure 3.30: Transient Response of Flexible Square Footing 


3.5.9 NATURAL FREQUENCIES OF A PARALLELEPIPED 


Code validation and mesh sensitivity studies for the BEST3D natural frequency capa- 
bility are being done using solutions obtained by Leissa and Zhang (Reference 33) for 
rectangular parallelpipeds. In the paper cited a Ritz technique is used to calculate the 
first five modes of each type (easy bending, stiff bending, torsion and extension) for a 
variety of geometries. Convergence studies in the paper indicate that the results for the 
first mode are accurate to better than 1%. 
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Figure 3.31: Transient Elasto- Plastic Response of a Bar Subjected to a Suddenly Applied 
and Maintained End Pressure 


One of the cases studied is a parallelepiped (Figure 3.32) of dimensions 1 x 1 x 0.5, with 
one end completely fixed. 

Four single region, quadratic BEST3D meshes (summarized in Table 3.5) were used in 
the study. The results of the study are shown in Table 3.6 for the first ten modes. Both 
the frequency parameter and the percentage deviation from the Ritz solution are given. 
It can be seen that rapid convergence is achieved for the first mode of each type by the 
addition of elements in the axial direction. Convergence of the higher modes, however, 
requires the use of additional elements in the chordwise direction as well. 

It is planned to use the NASA twisted plate data (Reference 34) and associated analyt- 
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Table 3.5: BEST3D Models for 1 X 1 X 0.5 Parallelepiped 


Model 

Elements 

Total Elements 

Source Points 


A 

B 

C 



i 

1 

1 

1 

6 

20 

2 

1 

2 

1 

10 

32 

3 

1 

3 

1 

14 

44 

4 

2 

3 

1 

22 
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Table 3.6: Convergence of BEST3D Natural Frequency Results 


Mode 

Type 

Ritz 

Mesh 1 

(%) 

Mesh 2 

(%) 

Mesh 3 

(%) 

Mesh 4 

(%) 

i 

EB 

0.447 

0.472 

5.5 

0.429 

-4.0 

0.435 

-2.7 

0.442 

-1.2 

2 

SB 

0.667 

0.664 

-0.4 

0.666 

-0.3 

0.668 

0.1 

0.661 

-0.8 

3 

T 

0.788 

0.887 

12.5 

0.829 

5.1 

0.820 

4.0 

0.788 

0.0 

4 

L 

1.596 

1.625 

1.8 

1.620 

1.5 

1.618 

1.4 

1.602 

0.4 

5 

EB 

1.664 

2.136 

28.3 

1.797 

8.0 

1.729 

3.9 

1.689 

1.5 

6 

SB 

1.774 



1.836 

3.5 

1.789 

0.9 

1.775 

0.1 

7 

T 

2.220 



2.552 

14.9 

2.448 

10.3 

2.285 

2.9 

8 

EB 

2.278 





3.033 

33.3 

2.365 

3.8 

9 

L 

2.797 







2.842 

1.6 

10 

SB 

3.068 







3.249 

5.9 

Mode Type Identification 


Frequency Parameter = 

Wyjp/E 


EB - 

Easy Bending 









SB - 

Stiff Bending 









T - Torsion 










L - Extension 










Table 3.7: Comparison of Ritz and BEST3D Results for Short, Thick Plate 


Mode 

Ritz 

BEST3D 

(%) 

1 (easy bending) 

3.38 

3.20 

-5.2 

2 (torsion) 

7.42 

7.20 

-2.8 
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Figure 3.32: Geometry for Calibration of BEST3D Natural Frequency Analysis 


ical results for the verification of the BEST3D natural frequency capability. The short, 
thick, untwisted plate (2 x 2 x .4) has already been analyzed using a model topologically 
equivalent to mesh 4 (above). The BEST3D results are compared, in Table 3.7, to results 
obtained using the Ritz method. The differences are somewhat larger than for the pre- 
vious case. This is believed to be due primarily to the degeneration in the worst element 
aspect ratios, from 1.5 for the earlier case to 2.5 for the NASA plate. This assumption 
will be verified by using a somewhat refined, two-region model. Analysis of the twisted 
plate geometries will then be undertaken. 


3.5.10 COMPONENT ANALYSIS 

1. Directionally Solidified Turbine Blade 

In order to verify the utility of the exact kernel functions (described in Ref. 2) for 
component analysis, the high turbine blade (Figure 3.33) previously analyzed with 
isotropic material properties was reanalyzed with the properties of a directionally 
solidified aerospace alloy. As expected, this static analysis ran without difficulty, 
indicating that both the kernel function for hexagonal materials and the new partic- 
ular solution for anisotropic centrifugal loads are operating correctly. The increase 
in cost relative to an isotropic analysis was approximately 30%. 

2. Bearing Race 

A single pitch of a bearing inner race was analyzed to determine the stress near a 
drain hole connecting the axial groove to the channel in which the bearing travels. 
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Figure 3.33: BEST3D Model of Anisotropic Turbine Blade 

The main load involved are a preload on the inner diameter, centrifugal load and 
the roller load. It is of interest to determine the response for a variety of loading 
combinations, and, in particular, for various locations of the roller relative to the 
axial groove. The BEST3D model used for this analysis is shown in Figure 3.34. It 
consists of three subregions containing a total of 323 surface patches. 

An initial analysis of the problem was carried out using linear variation of displace- 
ment and traction over all elements. As has been the case in previous complex 
analyses, the results were generally good but it was clear that improvement in 
accuracy was required in certain areas, particularly near the drain hole/groove in- 
tersection. To gain this improvement quadratic variation was used over 33 elements 
(10%). The increase in computing time was less than 25% and acceptable accuracy 
was achieved without the need to redefine the mesh. A typical contour plot of the 
hoop component of stress is shown in Figure 3.35. 
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Figure 3.34: BEST3D Model of Single Pitch of Bearing Race 





Figure 3.35: Hoop Stress Contours for Bearing Race 
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3.6 CONCLUSIONS 


A number of significant conclusion can be drawn from the experience of the first three 
years of this program. While the issues involved have been discussed in some detail above, 
the are presented here in summary form. 

1. The boundary element method is capable of producing accurate solutions to three- 
dimensional problems of plasticity, using a variety of algorithms. 

2. Preliminary results indicate that natural frequencies can be successfully calculated 
using the three-dimensional boundary element method. 

3. Elastodynamic analyses can be effectively treated using the boundary element 
method. While the nonlinear dynamic analysis has been formulated and imple- 
mented in preliminary form, a great deal of further work would be required to 
develop a usable tool for such problems. 

4. Reduction in volume discretization is a key to the efficient use of the boundary 
element method. Exploitation of particular integrals should be pursued in order to 
eliminate or reduce the need for volume cell definitions. 

5. Utilization of the special capabilities of supercomputers will be required for effective 
production use of three-dimensional boundary element analyses. 

6. Existing pre- and post-processing systems can be easily modified to interface with 
BEST3D. 
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3.7 LIST OF SYMBOLS 


Sii 

Kronecker delta symbol 

Ui,ti 

boundary displacements and tractions 

Gij 

displacement point load solution 

Fit 

fi 

traction kernel derived from Gij 
mechanical body forces 

T 

temperature 

(3 

coefficient of thermal expansion 

BijkyTijk , Dijki Sijk 

higher order kernels derived from Gij 

S 

surface of three-dimensional structure 

V 

interior of a three-dimensional structure 


stress tensor 

€ 'j 

strain tensor 

Cij 

jump terms in boundary integral equation 

X 

plastic flow factor 

Ffaj,h) 

yield function 

h 

hardening parameter 

*, same as & ij 

time derivative 

superscript F , as in 

ef . plastic component 

superscript e , as in 

ef. elastic component 

d ;,* 

elastic constitutive tensor 

x,y,z 

points in three-dimensional space 

superscript °, as in cr. initial stress (or strain) 


higher order kernel derived from Gij 

M 0 (v) 

isoparametric shape functions for volume cells 

* a (v) 

isoparametric shape functions for shape patches 

a b c 

oj ,u; , to 

Gaussian quadrature weights 

Va&c 

Gaussian quadrature points 

Cijki 

elastic constants for a general elastic material 

N,A 

6x6 matrices 


Lame constants 

t,T 

denote time in dynamic analysis 

*, as in Gij * U 

time convolution 

My, M 2 

shape functions for time variation 

[A^etc. 

coefficient matrices in time domain solution 
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